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0.1  Objectives 

The  objective  of  this  research  was  to  develop  tools  for  modeling,  analysis,  and  control  design  for 
unsteady  flow  and  combustion  phenomena  critical  to  operation  of  military  aeroengines. 


0.2  Summary  of  accomplishments 

We  have  made  significant  progress  in  control-oriented  modeling,  model  validation,  and  control  design 
for  processes  affecting  operation  of  aeroengines.  The  progress  in  modeling  includes  control-oriented 
reduced  order  models  for  turbomachinery  flow  including  compressors,  diffusers,  recirculation  zones, 
and  jets  in  cross  flow.  Major  breakthrough  was  achieved  in  the  area  of  nonlinear  model  validation 
and  parameter  identification,  with  particular  focus  on  models  of  combustion  instability.  In  particular, 
method  of  nonlinear  model  validation  and  parameter  identification  for  systems  with  non-equilibrium 
attractors  and  noise  drive  systems  based  on  comparison  of  long-term  statistics  has  been  transitioned 
to  internal  UTRC  projects.  New  results  in  control  of  limit-cycling  systems  with  bounded  control 
were  obtained.  Significant  progress  was  achieved  in  adaptive  control  of  combustion  instability,  control 
of  pressure  recovery  in  diffusers,  control  of  recirculation  zone,  and  control  of  jets  in  cross  flow.  In 
particular,  framework  for  control  of  mixing  based  on  hierarchy  of  vorticity  evolution  models  has  been 
introduced.  The  framework  has  been  used  in  control  of  mixing  enhancement  in  jets  in  cross-flow  and 
led  to  experimental  validation  of  model-based  analysis.  Results  of  the  research  are  in  13  conference 
papers  (10  published,  3  submitted)  and  7  journal  papers  (1  published,  4  submitted,  2  in  preparation). 
4  invited  sessions  were  organized  under  this  effort. 


0.3  Organization  of  the  report 

We  provide  extended  abstract  for  the  results  that  are  published  and  hence  publily  available.  We 
provide  more  details  for  the  results  that  are  not  available,  like  papers  that  are  submitted  or  in  the  case 
of  ungoing  research.  List  of  references  is  available  in  the  last  two  chapters  of  the  report.  References  to 
journal  papers  that  were  written  under  this  contract  are  referenced  with  letter  “j”  in  front  ([jl],  |j2], 
etc.).  Conference  papers  written  under  this  contract  are  referenced  with  letter  “c”  in  front  ([cl],  [c2], 
etc.). 
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Chapter  1 

Modeling  and  analysis  of  unsteady 
flows 

1.1  Adaptive  detection  of  instabilities  and  nonlinear  analysis  of  a 
reduced-order  model  for  flutter  and  rotating  stall  in  turboma¬ 
chinery 

In  paper  [cl]  by  Copeland,  Kevrekidis,  and  Ramiro  Rico-Martmez  system  parameters  were  adap¬ 
tively  varied  to  identify  bifurcations  from  equilibrium  in  simulations  of  a  reduced-order  model  for 
turbomachinery  aeromechanics.  An  element  of  the  adaptive  process  is  identification  of  a  low-order 
locally  nonlinear  discrete-time  map  from  the  observables.  The  nonlinear  behavior  of  the  system  in  the 
neighborhood  of  the  bifurcation  is  then  estimated  from  the  identified  system. 

Flutter  and  rotating  stall  are  instabilities  that  can  limit  the  design  and  operability  of  turbine 
engines.  Rotating  stall  is  primarily  aerodynamic  in  nature.  Flutter  involves  strong  coupling  between 
blade  vibration  and  unsteady  flow  perturbations.  Flutter  also  contributes  to  high-cycle  fatigue  and 
can  be  an  important  component  of  lifetime  cost.  In  contrast  to  rotating  stall,  the  current  capability 
for  flutter  prediction  is  poor,  hence  problems  are  often  found  late  in  development  cycle  and  remedial 
efforts  are  expensive.  Nonlinear  dynamics  may  be  important  because  rotating  stall  typically  occurs  as 
a  subcritical  bifurcation.  In  the  presence  of  subcritical  instability,  linear  analysis  alone  is  not  sufficient 
to  determine  boundaries  for  safe  operability. 

A  generic  approach  to  mapping  operability  regions  in  an  experiment  is  to  successively  set  the 
operating  parameters  to  fixed  values  and  passively  observe  the  system  dynamics  as  the  transients 
decay.  When  a  qualitative  change  in  the  long-term  behavior  of  the  system  is  observed  between  two 
consecutive  parameter  settings  it  may  be  inferred  that  a  bifurcation  exists  at  some  critical  intermediate 
parameter  value.  This  critical  value  may  be  bracketed  to  any  acceptable  accuracy  through  iterative 
bisection.  This  same  procedure  is  applied  generically  to  mathematical  models  when  the  only  tool  for 
investigation  is  direct  simulation.  This  is  often  the  case  in  computational  fluid  dynamics,  for  example. 

The  set-and-observe  approach  is  straightforward  but  it  has  some  obvious  shortcomings.  As  one 
approaches  a  bifurcation  point,  the  subcritical  state  of  the  system  becomes  less  stable  and  the  transients 
may  take  longer  and  longer  to  decay.  In  the  detection  of  hard  bifurcations  such  as  the  saddle-node 
bifurcation  or  the  subcritical  Hopf  bifurcation  there  is  the  additional  complication  that  overstepping 
the  critical  parameter  value  from  subcritical  to  supercritical  values  drives  the  system  far  from  the 
original  neighborhood  of  phase  space.  Bringing  the  system  back  to  the  neighborhood  of  the  bifurcation 
point  may  require  excessive  effort. 

In  this  paper,  detection  of  bifurcation  points  by  adaptively  varying  the  operating  parameters  is 
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demonstrated.  A  local  low-order  nonlinear  model  is  identified  from  input/output  data  as  the  system 
evolves.  This  identified  model  is  used  to  estimate  critical  parameter  values  for  bifurcation.  A  controller 
then  adjusts  the  parameters  (using  it  essentially  as  an  actuator)  to  bring  the  system  progressively  closer 
to  the  bifurcation  point.  The  identified  model  may  also  be  used  to  interpret  the  nonlinear  dynamics 
in  the  neighborhood  of  the  bifurcation. 

For  more  details  see  [cl]. 

1.2  Analysis  of  low  dimensional  dynamics  of  flow  separation 

The  increased  interest  in  the  modeling  and  active  control  of  separation  phenomena  [21,  107J  is  moti¬ 
vated  by  their  prevalence  in  several  industrial  situations  where  device  performance  is  often  hampered 
by  large  scale  separation  (e.g.  flow  within  engine  ducts  and  flow  past  high  angle-of-attack  airfoils). 
Until  recently,  separated  flows  have  been  studied  and  controlled  using  knowledge  about  the  hydro- 
dynamic  instabilities  that  govern  the  dynamics  near  the  separation  location.  These  studies  provide 
useful  insights  into  the  flow  physics  underlying  separation,  but  do  not  produce  models  that  can  be  used 
for  prediction  or  control.  The  use  of  the  dynamical  systems  approach  in  hydrodynamics  [4]  appears 
promising  for  the  development  of  reduced-order  models  of  flow  separation  dynamics.  In  particular, 
“parameter-dependent”  models  are  desired  that  can  be  used  to  predict  bifurcations  of  flow  states  and 
to  design  effective  control  methods. 

Low-dimensional  models  for  active  control  of  flow  separation  are  considered  in  the  papers  [c2]  by 
S.  Narayanan,  A.I.  Khibnik,  C.A.  Jacobson,  Y.  Kevrekedis  ,  R.  Rico-Martinez  and[jl]  by  K.  Lust  and 
A.I.  Khibnik,  S.  Narayanan,  C.A.  Jacobson,  and  K.  Lust. 

POD  method  was  used  to  analyze  the  dynamics  of  flow  separation  in  a  planar  diffuser.  The  pa¬ 
rameter  space  of  interest  is  associated  with  a  variable  diffuser  geometry,  given  here  by  the  expansion 
angle.  Two-dimensional  direct  numerical  simulations  of  the  diffuser  flow  were  performed  at  a  Reynolds 
number  of  5,200.  For  increasing  diffuser  angle,  the  flow  transitions  from  a  symmetric  steady  attached 
flow  to  an  asymmetric  separated  unsteady  flow  (also  termed  transitory  stall),  followed  by  a  steady 
separated  (stalled)  flow  and  leading  to  a  fully  separated  jet-like  flow.  Numerical  data  corresponding 
to  flow  fields  across  the  first  transition  (from  attached  flow  to  transitory  stall)  were  combined  to  yield 
a  set  of  POD  modes  that  span  the  dynamics  (both  transients  and  attractors)  for  different  expansion 
angles.  The  observed  dynamics,  inferred  from  the  first  few  POD  modes  and  their  coefficients,  included 
symmetry  breaking,  steady  states,  limit  cycles,  invariant  tori  and  chaotic  motion.  Limit  cycles  involve 
low  frequency  oscillations  associated  with  unsteady  stall  shedding,  independently  confirmed  via  labo¬ 
ratory  experiments.  A  bifurcation  scenario  (see  Figure  1.1)  was  hypothesized  to  provide  insight  into 
low  dimensional  dynamics  underlying  the  first  transition  and  to  guide  the  development  of  a  control 
strategy  for  separation.  The  POD  modes  were  also  used  to  construct  reduced  order  models  (Galerkin 
and  neural  network  models).  The  paper  briefly  discusses  these  approaches  and  a  proposed  discretized 
Galerkin  model. 

We  study  the  spatiotemporal  dynamics  of  flow  separation  in  a  planar  diffuser  to  extract  reduced- 
order  models  that  may  be  used  to  guide  active  separation  control.  Proper  orthogonal  decomposition  of 
numerical  simulation  data  revealed  organized  dynamics  of  the  large-scale  separated  structures.  This 
result,  combined  with  the  observation  (in  simulations  as  well  as  laboratory  experiments)  of  limit 
cycle/quasiperiodic  dynamics  and  spatial  symmetry  breaking,  suggests  an  underlying  low-dimensional 
dynamical  system.  Galerkin-based  and  neural  network-based  modeling  of  the  POD  modal  coefficients 
lead  to  (parameter  dependent)  low-dimensional  models  of  separation  dynamics.  The  use  of  such  models 
for  active  separation  control  appears  promising. 

Note  that  the  application  of  the  POD  for  the  modeling  of  deparated  diffuser  flows  has  not  been 
much  explored.  One  of  the  obstacles  that  makes  such  an  application  to  the  diffuser  non-trivial  is  the 
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Flow  classification  in  the  diffuser  based  on  flow  visualization  and  experimentally 
obtained  diffuser  pressure  recovery. 

The  corresponding  bifurcation  diagram  was  obtained  using  POD  modes  from  CFD 
computations. 


Bifurcation  scenario 


Aaycuictrte 
CMLwv  chvtic  ftjUx 


Figure  1.1: 


variable  geometry.  In  most  known  applications  of  POD  to  fluid  flows,  either  the  flow  conditions  or  the 
parameters  defining  the  boundary  conditions  are  assumed  to  be  fixed.  Computed  POD  modes  may 
still  be  used  to  span  the  flow  fields  for  different  flow  conditions  or  parameters  (e.g.,  varying  Reynolds 
number).  In  our  bifurcation  studies,  the  parameters  determine  the  geometry.  This  makes  it  difficult 
to  compare  flow  fields  at  different  angles  and  thus  to  compute  a  POD  basis  capable  of  approximating 
these  flow  fields.  To  address  this  issue,  we  transform  the  diffuser  geometry  to  the  “standard”  geometry 
of  a  channel  flow.  Thus,  the  parameters  enter  the  (transformed)  governing  equations  while  keeping 
the  geometry  and  the  boundary  conditions  fixed,  enabling  the  analysis  of  a  combined  set  of  flow  fields 
using  a  POD  basis. 

For  more  details,  please  refer  to  [c2jl]. 
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Chapter  2 
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* 


Nonlinear  model  validation  and 
parameter  identification 


2.1  Nonlinear  model  validation  and  parameter  identification  for  sys¬ 
tems  wifh  complex  dynamics  using  long-term  statistics  of  data 

2.1.1  Introduction 

In  [j2,c8]  we  presented  a  formalism  for  comparing  the  asymptotic  dynamics  of  dynamical  systems  with 
physical  systems  that  they  model.  There  is  often  no  need  for  the  detailed  (trajectory-wise)  comparison 
of  a  dynamical  system  and  the  physical  system  that  it  models,  but  only  comparison  in  statistical  sense. 
For  that  purpose,  invariant  measures  are  typically  considered.  But,  invariant  measures  usually  can 
not  be  observed  directly  in  an  experiment.  Thus,  we  base  our  formalism  on  time-averages  obtained 
from  a  single  observable  in  the  spirit  of  Takens  embedding  theorem.  In  particular,  we  constructively 
proved  that,  generically,  a  single  observable  is  needed  in  order  to  recover  an  invariant  ergodic  measure. 
Pseudometrics  on  space  of  dynamical  systems  can  be  defined  using  this  formalism  in  order  to  compare 
their  statistical  behavior.  We  also  identified  the  need  to  go  beyond  comparing  only  invariant  ergodic 
measures  of  systems  and  introduce  an  ergodic-theoretic  treatment  of  a  class  of  spectral  functionals 
that  allow  for  this.  We  discussed  the  extension  of  this  theory  to  stochastic  dynamical  systems. 

Both  the  concepts  of  invariant  measure  and  the  harmonic  average  formalism  developed  here  are 
related  to  spectral  properties  (in  particular,  the  point  spectrum)  of  the  so-called  Koopman  operator  U, 
a  linear  operator  that  acts  on  functions  on  the  phase-space  .  We  stress  that  in  this  context  questions 
of  identification  or  validation  of  asymptotic  properties  of  nonlinear  finite-dimensional  systems  with 
complex  dynamics  is  transferred  to  questions  of  identification  or  validation  of  a  linear,  albeit  infinite¬ 
dimensional  Koopman  operator.  Here  we  analyze  only  the  point  spectrum  of  that  operator.  Our  hope 
is  that  some  of  the  methods  developed  in  control  theory  of  linear  systems  can  be  used  to  study  these 
issues  further. 

The  ideas  introduced  in  [j2,c8]  can  be  used  for  parameter  identification  and  model  validation  of 
stochastically  driven  nonlinear  models  with  complicated  behavior. 

2.1.2  Example:  Parameter  identification  of  a  reduced-order  nonlinear  combustion 
model  with  noise  using  data  from  a  feedback  control  experiment 

As  an  illustration  we  provided  an  example  in  which  we  compared  the  asymptotic  behavior  of  a  combus¬ 
tion  system  measured  experimentally  with  the  asymptotical  behavior  of  the  model  that  is  a  stochastic 
control  dynamical  system. 
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Figure  2.1:  Control  phase  ramp  experiment.  Estimates  of  magnitude  and  frequency  of  pressure  oscil¬ 
lations  recorded.  Uncontrolled  mean  magnitude  level  indicated. 


Plane  2nd  order 
with  delay 


Valve  muration  Phase-shifting 
controller 


Figure  2.2:  A  model  of  a  controlled  combustion  process 

we  present  an  example  of  using  a  finite  number  of  functionals  evaluated  on  trajectories  to  identi¬ 
fication  of  parameters  of  a  model  describing  a  UTRC  combustion  rig  operating  close  to  an  unstable 
condition  associated  with  a  low  equivalence  ratio  condition.  In  an  experiment  a  feedback  loop  has  been 
closed  around  the  combustor  with  a  phase-shifting  controller.  The  control  phase  has  been  varied  and 
a  time  trace  of  pressure  has  been  recorded.  As  the  phase  varies,  both  attenuation  and  excitation  of  the 
pressure  oscillations  have  been  observed.  A  simple  model  of  the  controlled  combustion  process  consists 
of  a  Unear  plant,  controller,  a  saturated  actuator,  and  a  driving  broad-band  stochastic  disturbance. 
The  Unear  model  of  the  has  been  identified  from  an  experimentally  obtained  frequency  response.  The 
model  has  a  form  of  a  Ughtly  damped  2nd  order  system  with  delay.  An  analysis  of  the  linearization  of 
the  system  in  the  origin  indicates  that  the  Unear  eigenvalues  of  the  closed-loop  system  can  be  moved 
away  or  towards  the  imaginary  axis  on  the  complex  plane.  In  particular,  it  is  possible  to  move  the 
eigenvalues  to  the  right-half  plane.  In  this  case,  the  osciUations  would  grow  and  settle  on  a  Umit-cycle, 
as  the  actuator  saturation  prevents  an  unbounded  exponential  growth.  As  one  adds  the  strong  driving 
disturbance  in  the  model,  am  exact  analytical  study  of  the  behavior  of  the  system  is  no  longer  possible. 
One  can  stiU  run  a  numerical  simulation  and  compare  the  results  with  experiment.  We  assume  that 
the  disturbance  driving  the  combustion  model  is  a  white  Gaussian  process  with  standard  deviation  a. 
By  adjusting  the  standard  deviation  of  the  disturbance  and  the  value  of  delay  one  can  quaUtatively 
reproduce  the  results  of  the  phase  ramp  experiment  in  a  simulation.  In  what  foUows  the  delay  of 
the  plant  r  and  the  standard  deviation  of  the  disturbance  cr  wiU  be  varied  to  match  the  results  of 
experiment  with  closed  loop  control  for  several  values  of  control  phase.  We  will  define  a  measure  of  a 
good  fit  by  evaluating  several  long  term  statistics  and  combining  them  with  some  weights  into  a  one 
number. 

For  each  of  the  eight  values  of  the  control  phase  from  the  set  (-165,  -120,  -75,  -30,  15,  60,  105,  150) 


9 


Figure  2.3:  A  control  phase  ramp  obtained  from  a  model  simualtion. 


Figure  2.4:  Comparison  pf  PSDs  and  distributions  of  pressure  for  control  phase  9  =  -30.  Light  lines: 
from  experiment.  Dark  lines:  from  model  simulation.  Left:  for  r  =  0.006  and  a  —  1.  Right:  for 
r  —  0.006  and  o  =  3. 


(in  degrees)  we  chose  a  time  interval  containing  about  0.5  seconds  of  pressure  data  (1024  samples  of 
pressure  sampled  at  2000Hz)  corresponding  to  a  constant  control  phase. 

Below  we  show  the  time  traces,  phase  portraits,  distributions,  and  PSDs  of  pressure  for  control 
phases  9  =  —30  and  9  =  150  from  experiment  and  from  model  simulations  for  delay  r  =  0.006  and 
for  two  diffrent  values  of  the  standard  deviation  of  the  disturbance:  a  =  1  and  o  =  3.  The  control 
phases  6  =  -30  results  in  a  supression  of  the  pressure  oscillations.  The  corresponding  closed-loop 
model  is  a  stable,  disturbance  driven  system.  On  the  other  hand,  the  control  phase  6  =  150  results 
in  an  enhancement  of  oscillations.  The  model  analysis  without  disturbance  predicts  a  limit  cycling 
behavior. 

To  visualize  the  discrepancies  between  the  data  from  the  model  simulations  and  from  the  experi¬ 
ment,  we  present  the  PSDs  and  distributions  of  pressure  from  both  models  (darker  line)  on  the  same 
plot  with  the  experimental  PSDs  and  distributions  (lighter  lines).  We  can  see  that  the  lower  value 


.  i 


Figure  2.5:  Comparison  pf  PSDs  and  distributions  of  pressure  for  control  phase  9  =  150.  Light  lines: 
from  experiment.  Dark  lines:  from  model  simulation.  Left:  for  r  =  0.006  and  a  =  1.  Right:  for 
t  =  0.006  and  <7  =  3. 
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of  the  standard  deviation  of  the  disturbance  a  in  the  model  leads  to  an  underestimation  of  the  pres¬ 
sure  oscillation  level  in  experiment  for  both  control  phases.  We  also  observe  that  the  distribution  of 
pressure  in  the  data  from  the  model  simulation  for  high  value  of  standard  deviation  resembles  the  one 
from  experiment  better  than  the  data  from  the  low  disturbance  model  simulation. 

To  quantify  the  discrepancies  of  the  data  from  experiment  and  model  simulations  we  would  like 
to  choose  some  quantities  that  describe  the  average  frequency  of  the  oscillations,  the  average  size  of 
the  oscillations,  and  somehow  distinguish  a  disturbance-driven  limit  cycle  from  a  disturbance-driven 
stable  system.  With  this  goal  in  mind,  for  each  time  interval  corresponding  to  a  given  control  phase 
6  we  calculated  the  following  measures  of  the  long  term  statistics: 

1.  The  mean  frequency  of  the  pressure  oscillations. 

2.  The  standard  deviation  of  the  frequency  pressure  oscillations. 

3.  The  standard  s  deviation  of  the  pressure  signal. 

4.  The  ratio  of  pressure  measurements  in  the  interval  [-s/2  s/2]  to  the  number  of  all  pressure 
measurement  in  a  given  time  interval. 

(In  the  sequel  we  will  refer  to  these  measures  as  the  long  term  statistics ). 

The  mean  frequency  and  the  standard  deviation  of  the  frequency  of  pressure  oscillations  has  been 
calculated  from  the  Power  Spectral  Density  of  the  pressure  signal.  The  mean  frequency  of  the  pressure 
oscillations  and  the  standard  deviation  of  the  pressure  signal  are  of  obvious  interest.  For  a  stable  linear 
system  the  standard  deviation  of  the  frequency  is  a  measure  of  the  damping  of  the  system.  The  ratio 
of  pressure  measurements  in  the  interval  [-s/2  s/2]  to  the  number  of  all  pressure  measurement  in 
a  given  time  interval  was  chosen  to  distinguish  between  a  disturbance-driven  stable  system  from  a 
disturbance-driven  linearly  unstable,  limit-cycling  system.  As  we  will  see,  on  the  stable  side  this  value 
is  typically  around  .35-. 4,  while  on  the  unstable  side  it  drops  below  .3. 

A  measure  of  a  discrepancy  of  the  statistics  from  experimental  data  and  model  simulation  was 
obtained  by  taking  the  absolute  value  of  the  difference  of  the  statistics  and  dividing  it  by  the  value 
of  the  statistic  for  the  experimental  data.  A  mean  value  of  this  relative  measure  for  all  eight  control 
phases  is  called  in  the  sequel  a  relative  error  for  a  given  statistic. 

From  the  four  relative  errors  we  construct  a  total  error  measure  by  adding  the  four  relative  errors 
with  the  following  weights: 

1.  The  weight  for  the  mean  frequency  of  the  pressure  oscillations.  0.3 

2.  The  weight  for  the  standard  deviation  of  the  frequency  pressure  oscillations.  0.1 

3.  The  weight  for  t^e  standard  s  deviation  of  the  pressure  signal.  0.4 

4.  The  weight  for  the  ratio  of  pressure  measurements  in  the  interval  [—s/2  s/2]  to  the  number  of  all 
pressure  measurement  in  a  given  time  interval.  .2 

The  values  corresponding  to  the  lowest  total  relative  error  of  0.09321  are  r  =  0.006  and  a  =  3.  The 
comparison  of  the  long  term  statistics  from  experiment  and  model  simulation  for  this  case  is  shown 
in  Figure  2.6.  Note  that  the  long  term  statistics  introduced  in  this  section  to  select  a  combustion 
model  that  best  matches  data  were  chosen  in  an  arbitrary  manner,  to  measure  various  features  that 
were  considered  of  merit.  The  results  of  this  paper  were  not  a  necessary  for  being  able  to  pick  certain 
statistics  in  this  particular  example.  What  we  would  like  to  stress  here  is  that  all  these  different  and 
somewhow  unrelated  statistics  are  constructed  as  certain  functionals  defined  on  recorded  trajectories 
from  an  experiment  and  from  a  model  simuation.  The  results  of  this  paper  provide  a  new  insight  into 
what  is  the  precise  mathematical  meaning  of  often  ad  hoc  chosen  measures  of  a  distance  between  two 
dynamical  systems. 

A  single  peak  in  PDF  of  pressure  indicates  a  stable  driven  system.  Double  peak  indicates  a  noisy 
limit  cycle.  An  alternative  method  of  distinguishing  whether  a  model  that  is  a  noise-driven  stable 
system  or  a  noise-driven  limit  cycling  system  better  mathces  experimental  data  is  presented  in  [36]. 
The  method  proposed  in  [j2,c8]  became  a  standard  method  for  validation  of  reduced  order  models 
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Figure  2.6:  Long  term  statistics  from  experiment  and  model  simulation  for  r  =  0.006,  a  =  3. 

of  combustion  instability  used  at  UTRC.  In  particular,  the  method  was  used  in  an  internally  funded 
combustion  programs  for  model  validation  and  parameter  identification  a  single-nozzle  rig  model. 
Figure  2.7  shows  comparison  of  Power  Spectral  Density  and  Probability  Density  Function  of  the 
pressure  data  from  experiment  and  model  simulation.  Moreover,  the  method  has  been  applied  to 
identification  of  a  reduced  order  model  of  combustion  instability  in  an  advanced  military  engine  and 
to  distinguish  regions  of  validity  of  stable  driven  and  limit  cycling  models  of  military  aeroengines  and 
industrial  gas  turbines.  The  method  can  be  also  used  for  validation  of  reduced  order  models  of  shear 
flows. 

2.2  Nonlinear  model  validation  using  frequency  reponse  data 

In  [c4]  we  presented  a  model  validation  problem  using  a  frequency  response  data,  motivated  by  the 
modeling  of  a  combustion  process  at  UTRC.  The  system  excitation  was  well  modeled  by  a  stochastic 
source  and  this  dominated  the  response.  The  number  of  data  points  that  must  be  recorded  to  capture 
the  system  behavior  precluded  the  use  of  tune-domain  model  validation  techniques.  A  critical  analysis 
of  the  results  of  the  linear  model  validation  indicated  the  limitations  of  the  linear  model  validation 
techniques  in  the  context  of  nonlinear  feedback  systems.  The  need  for  new  model  validation  tech¬ 
niques  that  take  into  account  nonlinear  model  perturbations  and  the  effect  of  the  driving  broad  band 
disturbance  was  indicated. 

Frequency  response  measurements  of  UTRC  single  nozzle  rig  combustor  were  taken  in  both  open- 
loop  and  closed-loop  configurations  using  swept-sine  experiments.  The  results  of  these  experiments 
are  shown  in  Figure  2.8.  In  both  cases  a  very  good  fit  to  the  response  was  obtained  with  a  stable 
linear  transfer  function  consisting  of  a  2nd  order  system  with  delay. 

The  nonlinearity  assumed  in  the  model  did  not  manifest  itself  by  distorting  the  experimental 
frequency  response  estimates.  As  a  consequence,  it  is  natural  to  quantitatively  investigate  the  extent 
to  which  a  linear  model — with  a  norm  bounded  perturbation — can  describe  the  system  behaviors. 
Such  a  model  structure  is  illustrated  in  Figure  2.9. 

Both  linear  and  nonlinear  models  contain  bounded  perturbations,  which  allow  each  to  describe 
a  set  of  behaviors.  This  structure  is  commonly  used  in  robust  control  and  provides  a  degree  of 
robustness  with  respect  to  the  designer’s  uncertainty  about  the  true  range  of  system  behaviors.  Given 
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Figure  2.7:  Results  of  model  parameters  fitting  using  statistics  of  data  from  experiment  and  model 
simulation. 
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£££ , ,^pP"  &eqUenCy  0t  ^  — d  *°  T”  under  open-bop 


experimental  input-output  data,  model  validation  techniques  provide  a  means  of  comparing  the  ability 
of  each  model  to  account  for  the  observed  experimental  data. 

closed-loop  response  appeared  linear,  it  was  significantly  different  from  that 

th  t6  1°“  6  °pen-loop  modeL  ^  the  closed-loop  response  the  resonant  peak 

"  ^  the  °Pen-lo°P  case,  and  the  damping  was  half  that  of  the  open-loop  case  (refer  to 

nffhp  a'8  ‘  iTh  >  rue  Tr  h0Ugh  the  pressure  oscillations  in  the  closed-loop  case  were  half  those 
of  the  open-loop  situation.  A  hnear  analysis  of  the  control  design  indicates  that  we  would  expect  the 
damping  to  increase  under  closed-loop  control.  V 

An  interpretation  of  this  unexpected  behavior  can  be  made  based  on  the  nonlinear  model  as  follows 
A  strong  drrnng  noise  adds  to  the  sinusoidal  driving  signal  at  the  input  of  the  nonlinear  function  and  a 

Xnr  f!rTtte?. part  the  fUnCtiOD  “  Visited-  111  a  closed-loop  experiment  the  controUer 
reduces  the  size  of  the  disturbance  and  less  of  the  saturated  portion  of  the  nonlinear  function  is  visited. 

Effectwely,  the  nonlinearity  provides  a  higher  gain  to  the  sinusoidal  signal  and  hence  the  positive 
feedback  m  he  internal  loop  is  stronger.  This  manifests  itself  as  a  higher  resonant  peak  and  lower 
damping  m  the  hnear  model  fit  to  the  combustor  frequency  response  in  the  closed-loop  experiment. 
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Figure  2.9:  Linear  perturbation  model  of  the  combustion  process 


Figure  2.10  illustrates  a  model  structure  capable  of  reproducing  the  observed  behavior. 

The  closed-loop  system,  illustrated  in  Figure  2.10,  has  been  simulated  and  used  to  obtain  the 
simulated  frequency  response  measurements  shown  in  Figure  2.11.  The  simulation  is  a  very  good 
qualitative  match  to  the  observed  experimental  behavior:  the  damping  ratio  of  the  oscillatory  mode 
is  reduced  in  both  the  closed-loop  model  simulation  and  experiments. 


2.3  Nonlinear  model  identification  using  frequency  reponse  data 

2.3.1  Introduction  and  problem  formulation 

In  papers  [j7,cl3],  we  propose  a  describing  function  based  procedure  for  identifying  the  shape  of  the 
nonlinearity  in  the  thermoacoustic  feedback  loop.  We  describe  the  procedure  in  two  cases: 

1.  where  no  driving  noise  is  present  and  the  system  is  either  stable  or  in  self-sustained  limit  cycling 
oscillations  and 

2.  where  Gaussian  driving  noise  alone  is  present  and  limit  cycle  is  absent. 

In  the  case  where  no  driving  noise  is  present,  we  show  that  the  nonlinearity  shape  may  be  recovered 
by  using  experiments  with  feedback  control  where  the  control  is  used  to  vary  the  amplitude  of  the 
limit  cycle  from  zero  to  the  maximum  value.  At  each  fixed  amplitude,  standard  linear  identification 
is  used  for  identifying  the  harmonic  balance  gain.  These  gains  are  then  used  to  determine  tl  i  shape 
of  the  nonlinearity  function. 

In  the  case  where  Gaussian  driving  noise  is  present,  a  procedure  is  proposed  whereby  an  external 
controller  is  used  to  manipulate  the  variance  of  the  Gaussian  noise  seen  by  the  nonlinearity.  At  each 
controller  setting,  small  amplitude  sinusoidal  forcing  is  used  to  determine  the  harmonic  balance  gain 
(in  this  case,  for  the  Gaussian  noise).  The  shape  of  the  underlying  nonlinearity  is  recovered  from  these 
harmonic  balance  gains. 

A  simplified  version  of  UTRC  combustion  model  of  a  premixed  combustion  process  is  an  intercon¬ 
nection  of  a  linear  lightly  damped  oscillator,  delay,  and  a  saturating  nonlinearity  in  a  feedback  loop. 
Figure  2.12  gives  the  block  diagram  schematic  of  the  combustion  model.  The  oscillator  represents  the 
acoustic  waves  in  the  combustion  chamber.  The  delay  represents  the  time  it  takes  the  fuel  and  air 
to  mix  in  a  premixing  nozzle,  transport  to  the  flame  front,  and  react.  The  saturating  nonlinearity 
represents  the  dependence  of  a  heat  release  rate  on  the  fuel  to  air  ratio  at  the  flame  front.  The  feed¬ 
back  loop  represents  the  effect  of  oscillations  of  air  mass  flow  in  the  nozzle  on  the  fuel  to  air  ratio  and 
hence  on  the  heat  release  rate.  The  delay  is  such  that  the  internal  feedback  loop  provides  a  positive 
feedback  around  the  oscillator.  Only  the  output  of  the  oscillator  (pressure)  is  accessible  for  measure¬ 
ment.  Depending  upon  the  experimental  conditions,  a  strong  broad  band  disturbance  (representing 
turbulent  velocity  fluctuations)  may  be  driving  the  system  at  the  input  of  the  linear  oscillator.  The 
control  input  (representing  fuel  modulation)  adds  to  the  disturbance  input. 
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Figure  2.10:  Nonlinear  model  structure  matching  the  observed  closed-loop  experimental  frequency 
responses. 
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Figure  2.12:  Nonlinear  model’s  block  diagram 

The  nonlinear  identification  objective  is  to  obtain  the  shape  of  the  heat  release  nonlinearity  by  using 
the  control  and/or  the  external  forcing  (see  Figure  2.12).  For  the  purpose  of  this  note,  it  is  assumed 
that  the  linear  dynamics  (acoustics  and  delay)  have  already  been  obtained  using  a  frequency  response 
test  and  are  known.  Additionally,  we  impose  a  special  structure  on  the  heat  release  nonlinearity  as 
described  by  the  following  definition. 

Definition.  The  heat  release  nonlinearity  h{u)  is  a  static  function  of  its  input  u  such  that  the  following 
three  requirements  are  satisfied 

1.  h{0)  =0, 

2.  h  is  a  non-decreasing  odd  function  ( h(—u )  =  —h(u))  and 

3.  For  u  >  0,  h(u)  is  a  concave  function,  i.e. 

h(Xu  +  (1  -  AM  >  A/i(u)  +  (1  -  X)h(v).  (2.1) 

For  such  a  non-linearity  and  for  u  >  0,  we  have 

h(Xu)  >  A/i(u),  0  <  A  <  1  (2.2) 

A h(u)  >  h{ Xu),  A  >  1.  (2.3) 

We  will  have  the  need  to  use  these  inequalities  in  certain  proofs  to  follow.  A  relay  nonlinearity  is  an 
example  of  heat  release  nonlinearity. 

Remark.  If  additionally,  h  is  assumed  to  be  a  C1  function  then  the  requirement  (3)  in  the  definition  of 
the  heat  release  nonlinearity  can  be  shown  to  be  equivalent  to  the  derivative  h!  being  a  non-increasing 
function  of  its  argument. 
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2.3.2  Identification  without  noise 

In  the  absence  of  noise,  harmonic  balance  for  the  limit  cycling  system  yields  the  sinusoidal  describing 
function  gain 

1  ft* 

NA  =  -h  h{Asin{0))sin{e)d6  (2.4) 

IT  A  Jo 

for  the  static  nonlinearity  h  where  A  represents  the  amplitude  of  the  sinusoid  at  the  input  of  the  heat 
release  nonlinearity.  This  sinusoid  may  arise  either  due  to  the  external  forcing  or  due  to  the  limit 
cycling  oscillations.  In  the  limit  cycling  case,  the  loop  equation 

1  +  NAG0(iu))e-(iwr)  =  0  (2.5) 

needs  to  be  additionally  satisfied.  In  this  section,  we  do  not  consider  the  case  where  both  the  forcing 
and  limit  cycle  oscillations  are  present. 

If  the  thermoacoustic  system  is  non  limit  cycling  and  the  presence  of  the  sinusoidal  forcing  does  not 
cause  the  system  to  exhibit  limit  cycling  behavior  then  harmonic  balance  may  be  used  together  with  the 
sinusoidal  forcing  to  estimate  the  feedback  gain  due  to  the  nonlinearity  (at  a  fixed  forcing  amplitude). 
By  estimating  this  gain  for  many  forcing  amplitudes,  the  harmonic  balance  gain  is  computed  for 
different  amplitudes  at  the  nonlinearities  input. 

If  the  thermoacoustic  system  is  limit  cycling  then  the  same  procedure  is  applied  but  this  time  with 
an  external  controller.  The  controller  is  gain  and  phase  adjusted  to  time  the  amplitude  of  the  limit 
cycle.  The  harmonic  balance  gain  is  now  evaluated  at  the  different  limit  cycle  amplitudes. 

To  compute  the  nonlinearity  shape  from  these  gains  identified  using  linear  techniques,  we  arrive 
at  the  operator  equation 

1  r2* 

g(A)  =  Lh(A)  =  —  h(Asin(6))sin(6)dO  (2.6) 

7T  Jq 

which  needs  to  be  solved  to  determine  the  shape  of  the  nonlinearity.  For  the  heat  release  nonlinearity 
(see  definition  in  section  2.3.1),  this  equation  is  invertible  and  the  inverse  is  given  by 

h(A)  =  '£Mng(A)  (2.7) 

n=Q 

where  the  operator 

M  =  I  —  L.  (2.8) 

The  proof  that  the  series  actually  converges  for  the  heat  release  nonlinearity  requires  ||M||  to  be  less 
than  1.  This  follows  because 


Mh(A) 


(I  —  L)h(A) 
h{A)  -  Lh{A) 

—  f  h(A)sin2(9)d6  —  —  [  h(Asin(6))sin(0)d6 

Jo  a-  Jo 

—  [  [h(.A)stn(0)  —  h{Asin(9))]  sin(6)d9 

7T  Jo 


(2.9) 


because  h  is  assumed  to  be  an  odd  function.  Finally,  because  h  is  concave,  we  use  the  inequality  (2.2) 
to  estimate  a  bound  for  the  operator  as 


\Mh(A)\  < 


-  I**  [1  -  sin(0)]  h{A)sin(e)de 
TT  Jo 

(;  -  DAM)- 

7 r 


(2.10) 


19 


(2.11) 


We  thus  have  a  bound  on  the  operator  norm  of  M  as 

1WII  <(;  —  »)<  1- 

7 r 

A  Matlab  code  has  been  written  to  compute  the  terms  in  the  series  iteratively  (see  Appendix).  The 
series  representation  of  the  inverse  also  allows  one  to  obtain  bounds  on  the  error.  For  example,  a  relay 
type  nonlinearity 


h{A)  =  6,  A  >  0 
=  -6,  A  <  0 

yields  through  a  harmonic  gain  analysis  Na  —  or  equivalently 

r2ir 


1  rZlT  4  b 

glA)  =  -  /  bsin{9)de  =  — . 

TT  JO  K 


(2.12) 


(2.13) 


The  computations  are  then  easily  performed  to  recover  the  shape  of  the  nonlinearity.  The  0(0) 
computation  yields 

A\  = 


4  b 


with  relative  error  bounded  by 


h°(A)  =  g(A)  =  -  =  1.27236 

7 r 


<  IWI  <  (|  - 1). 

\\h\\  7T 


The  0(1)  computations  yield 

h\A)  =  g{A)  +  Mg(A) 


46  2  n  [46  .  ...  461  . 

= - 1 —  /  — stn(9) - stn(9)d9 

7T  n  Jo  IT  T  . 


=  0.92536 


with  relative  error  bounded  by 


Tf#  <  IIMII2  <  (i  -  l)2  =  0.0747. 

I \h\\  * 


Similarly  relative  error  for  0(2)  compuations  is  bounded  by  0.0204  and  so  on. 


(2.14) 

(2.15) 


(2.16) 


(2.17) 


2.3.3  Identification  with  noise 

In  the  presence  of  Gaussian  noise,  the  random  input  describing  function  framework  is  used  to  do 
harmonic  balance.  Throughout  this  section,  we  assume  that  limit  cycle  is  not  present  and  the  harmonic 
balance  due  to  the  noise  alone  accounts  for  the  signals  in  the  loop.  If  the  variance  of  the  signals  at 
the  input  and  the  output  of  the  nonlinearity  is  known,  the  operator  equation 

g(a)  —  Lh{a)  =  —==.  f  h(cs)se~s  ! 2ds  (2.18) 

V  27T  J — oo 

needs  to  be  solved  to  obtain  the  nonlinearity  shape  h.  For  the  heat  release  nonlinearity,  the  solution 
may  again  be  obtained  in  a  series  form 

OO 

Mct)  =  Mng(&)  (2-19) 

n=0 
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where  the  operator 


M  =  I  —  L.  (2.20) 

For  the  heat  release  nonlinearity  (see  Definition  in  section  2.3.1)),  the  inequalities  (2.2)  and  (2.3)  may 
be  used  to  carry  out  an  analysis  similar  to  the  one  described  in  previous  section  to  show  that 


\Mh(c)\  =  \\r~  !  [fi(<r)s  -  fi(crs)]  se  s2/2ds| 

V  7T  Jo 

<  | \J~  J  [h{cr)s-h(cTs)]se~s,/2ds\  +  \}J^  [h(<r)s  -  h(as)]  se~si^ds\ 

<  (/  (1  “  s)se~*2/2ds  +  (s  -  l)se-a2/2ds^ 

=  0.435||^|| 

<  l||fc||. 

For  the  relay  type  nonlinearity,  the  harmonic  balance  gain  Nr  —  £  and  equivalently 

9 (<r)  =  ~J==  f  bse~3i/2ds  =  ^6. 

v27T  J  —  oo  Tf 

The  series  evaluation  for  the  inverse  then  gives 

h°(a)  =  g(a)  =  —  b  ss  0.86, 

7 r 

/i*( a )  =  g(a)  +  Mg{a)  «  0.966. 


(2-21) 


(2.22) 


(2.23) 


(2.24) 

(2.25) 


with  corresponding  relative  errors  bounded  by  0.435  and  0.189  respectively. 

In  combustion  applications  the  variance  of  the  signals  both  at  the  input  and  the  output  of  the 
nonlinearity  is  not  known  (as  heat  release  rates  can  not  typically  be  directly  measured).  As  a  result, 
the  procedure  described  above  may  not  be  directly  applied.  In  the  following  paragraphs,  we  describe 
a  procedure  for  estimating  the  heat  release  nonlinearities  random  input  describing  function  gain  by 
using  small  magnitude  sinusoidal  forcing.  Before  describing  the  procedure,  we  state  a  key  result  which 
makes  the  estimate  feasible. 

Result  For  a  heat  rdaase  nonlinearity,  the  random  input  describing  function  gain 

NR(*)=KmNA(cT).  (2.26) 

y4->0 

If  the  function  h  is  additionally  C1,  this  gain  is  given  by 

Nr(c)  =  -L  f°°  h'{<7 s)e~(,2/2)ds.  (2.27) 

v27T  J  —  oo 

Proof.  We  first  assume  the  function  h  to  be  C1.  The  random  input  describing  function  gain  (in  the 
absence  of  any  limit  cycle  oscillations)  is  given  by 


Nr(ct,A  =  0)  =  --L=  f  ^^-se  *3/2ds. 

V  27T  J  —  oo  g 


y/2 n 

Using  integration  by  parts,  we  may  express  the  right  hand  side  as 


(2.28) 


Nr 


1  h(os)__t,  2 
"  V2t  c 

=  1 —  [  h'(crs)e~s7/2ds. 

V2n  7-oo 


00  1  ro o 
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(2.29) 
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This  gain  needs  to  be  estimated  using  sinsuoidal  forcing.  The  sinsusoidal  describing  function  gain  (in 
the  presence  of  noise  of  variance  a) 

Na{c,  A)  =  /l  f  d9  [  dsh{Asin{9)  +  as)sin(9)e~s2/2.  (2.30) 

(27T)'3' Jo  j-oo 


Since  Jq  sin(0)d6  -  0,  we  may  write  this  gain  as 


NA(c,A)  =  *  J\e  jyV±l±^ftL USsAUmc**.  (2-31) 


(2tt)3/2 

Now  applying  Dominated  Convergence  Theorem  we  have 


O  f  27T  f  oo  2 . 

lim  Na(<t,A)  =  7—7377  d0  ds  ti{as)sin2{9)e~5  1 
A->0  (27r)3/2  Jo  J-OO 

=  -J=  f°°  h '{as)e~s2'2ds, 

v/2W-oo 


(2.32) 


The  equations  (2.29)  and  (2.32)  together  give  the  result  for  a  C 1  nonlinearity.  For  arbitrary  bounded 
nonlinearity  function  h,  we  simply  approximate  the  nonlinearity  by  a  sequence  of  C1  functkus  which 
converge  (in  a  pointwise  sense)  to  the  function  h  and  apply  DCT  to  obtain  the  result. 

The  significance  of  the  result  is  that  the  random  input  describing  function  gain  may  be  estimated 
by  carrying  out  a  frequency  response  test  with  vanishing  sinusoidal  signal  amplitude.  Of  course,  for  a 
system  forced  with  very  large  bandwidth  noise,  one  would  get  in  to  problems  related  to  SNR  so  any 
estimate  of  the  noise  gain  Nr  would  necessarily  be  biased. 

Remark.  The  exact  expression  for  Nr,  the  random  input  describing  function  gain  also  shows  that  as 
the  variance  c  of  the  Gaussian  noise  at  the  input  of  the  heat  release  nonlinearity  decreases,  the  gain 
Nr(cf)  increases. 

The  random  input  describing  function  gain  is  estimated  by  using  a  frequency  response  test.  A 
feedback  controller  is  used  to  manipulate  the  variance  of  the  Gaussian  noise  seen  at  the  input  of  the 
heat  release  nonlinearity.  For  each  setting  of  the  controller,  a  sinusoidal  input  is  used  to  estimate 
the  random  input  describing  function  gain  Nr  and  this  together  with  the  variance  computed  from 
the  time  series  is  used  for  estimating  the  variance  of  the  signals  at  the  input  and  the  output  of  the 
nonlinearity.  By  procuring  these  variances  at  different  control  settings,  sufficient  data  is  obtained  to 
apply  the  procedure  to  compute  the  nonlinearity  shape. 


2.3.4  Simulations 

A  simulink  model  was  used  to  build  a  mathematical  model  for  the  combustion  block  diagram  described 
in  Figure  2.12.  We  describe  two  simulations  which  were  carried  out  to  demonstrate  the  procedure 
described  in  sections  2.3.2  and  2.3.3. 

First,  we  consider  the  no  noise  case.  In  the  absence  of  noise,  the  combustion  is  unstable  and  a 
limit  cycling  system.  The  controller  is  used  to  vary  the  linear  dynamics  of  the  loop  and  the  resulting 
amplitude  of  the  pressure  time  series  oscillations  and  the  control  input  time  series  oscillations  observed. 
These  two  amplitudes  together  with  their  phase  differences  and  information  on  the  linear  dynamics 
present  in  the  loop  are  used  to  estimate  the  amplitude  information  at  the  input  and  the  output  of  the 
nonlinearity.  These  estimates  are  then  used  to  estimate  the  nonlinearity  shape.  Figure  2.13  compares 
the  actual  nonlinearity  with  the  one  obtained  from  the  experiment. 

Next,  we  consider  the  case  where  a  strong  Gaussian  noise  component  is  present.  In  the  absence 
of  control,  Gaussian  noise  balance  accounts  for  the  signals  in  the  thermoacoustic  loop  (actually,  a 
small  limit  cycle  is  present  -  see  Figure  2.15).  The  control  is  used  to  manipulate  the  variance  of  the 
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Nonlinearity  shape  x  (exact)  and  o  (estimate) 


Figure  2.13:  Estimate  of  the  nonlinearity  from  the  harmonic  gains  -  limit  cycling  case 

signal  as  seen  at  the  input  to  the  nonlinearity.  A  small  amplitude  forcing  may  be  used  to  estimate 
the  nonlinearity  gain  at  each  control  setting  as  described  before.  The  objective  of  this  forcing  test  is 
to  estimate  the  random  input  describing  function  gain.  However,  we  have  not  been  able  to  succesfully 
carry  out  this  portion  of  the  simulation  yet. 

Instead,  we  use  the  signals  at  the  input  and  the  output  of  the  nonlinearity  directly  to  estimate  the 
harmonic  balance  gains.  Figure  2.14  compares  the  nonlinearity  shape  estimated  with  this  procedure 
versus  the  actual  nonlinearity.  The  error  between  the  two  is  due  to  the  fact  that  the  input  and  the 
output  signals  are  not  pure  Gaussian  signals  rather  are  a  combination  of  Gaussian  and  a  sinusoidal 
limit  cycle.  Figure  2.15  plots  the  PDF  (histograms)  for  the  pressure  time  series.  A  single  peak  in  PDF 
of  pressure  indicates  a  stable  driven  system.  Double  peak  indicates  a  noisy  limit  cycle.  For  stabilizing 
control  action  (positive  control  gains),  the  signals  are  Gaussian  and  the  PDF  does  not  display  any 
limit  cycle  component  while  for  destabilizing  gains(negative  control  gains),  the  limit  cycle  becomes 
prominent. 

If  one  were  to  blindly  treat  the  signals  as  Gaussian  (as  we  have  done  here),  the  variance  estimate  is 
biased  and  the  harmonic  balance  gains  thus  derived  are  incorrect.  By  the  nature  of  the  procedure,  the 
error  is  larger  for  large  a  where  the  control  is  most  destabilizing  (and  hence  the  limit  cycle  component 
is  the  largest).  A  partial  fix  to  this  problem  is  to  use  Prabhir’s  algorithm  to  isolate  the  a  for  the 
input  and  the  output  signals  and  use  these  instead.  This  will  provide  an  improved  estimate  of  the 
nonlinearity  shape  if  the  limit  cycle  amplitude  is  not  significant. 

However,  a  complete  solution  of  the  problem  where  both  the  limit  cycle  and  the  Gaussian  noise 
components  are  together  present  requires  an  identification  algorithm  which  incorporates  the  multiple 
input  describing  function.  Though,  this  will  be  subject  of  future  work,  next  section  discusses  some  of 
the  issues. 

2.3.5  Limitations  of  the  identification  procedure 

The  available  control  authority  fundamentally  limits  the  identification  of  the  nonlinearity.  In  the 
simulations  provided  here,  unlimited  control  authority  was  assumed.  However,  the  effect  of  limited 
control  authority  for  the  Gaussian  noise  and  the  limit  cycling  cases  are  fundamentally  different. 

For  the  limit  cycling  case,  if  one  has  information  on  harmonic  balance  gains  for  limit  cycle  ampli- 


23 


12 


Nonlinearity  shape  x  (exact)  and  o  (estimate) 


Figure  2.14:  Estimate  of  the  nonlinearity  from  the  harmonic  gains  -  Gaussian  noise  case 
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Figure  2.15:  PDF  of  the  pressure  time  series  (input  to  nonlinearity)  for  different  control  gains  -  positive 
control  gains  are  stabilizing  and  negative  destabilizing 
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tudes  from  0  to  A  then  one  can  (in  theory)  identify  the  nonlinearity  shape  for  all  values  of  input  from 
0  to  A  with  out  any  error.  This  is  because  of  the  nature  of  the  operator  M  in  the  case  of  limit  cycling 
(or  sinusoidally  forced)  systems. 

For  the  noise  driven  system  however,  in  order  to  identify  the  nonlinearity  precisely  at  any  input 
level,  information  on  harmonic  balance  gains  from  0  to  oo  is  required.  This  information  is  impossible 
to  obtain  from  any  realistic  experiment.  In  any  practical  implementation  of  the  procedure,  the  data 
points  outside  the  range  available  are  interpolated  to  recover  the  shape  of  the  nonlinearity.  As  a  result, 
any  identification  with  Gaussian  noise  is  necessarily  prone  to  errors. 

Finally,  we  discuss  the  problem  where  both  the  limit  cycles  and  the  Gaussian  noise  are  together 
present.  The  problem  there  is  to  back  out  the  nonlinearity  shape  from  informations  on  Nr  and  Na 
(the  harmonic  balance  gains  for  the  noise  and  limit  cycles  respectively)  together  with  assoscb.ted  a 

and  A. 
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Chapter  3 


Control  of  non-equilibrium  behavior  in 
aeroengines 


3.1  Control  of  planar  limit  cycles 

This  section  is  an  extended  abstract  of  results  presented  in  paper  [c3]. 

There  are  numerous  examples  of  systems  where  oscillations  are  attributed  to  stable  limit  cycles. 
Reduced  order  models  of  combustion  instability  [93]  and  compression  system  dynamics  [37]  are  in  this 
category.  The  control  goal  is  to  reduce  the  size  or  elimination  of  the  limit  cycle.  One  group  of  methods 
for  control  design  are  linear  methods  for  damping  augmentation  of  the  unstable  origin.  The  methods 
axe  essentially  the  same  as  the  ones  for  damping  augmentation  of  the  stable  linear  systems.  With 
large  control  authority  the  origin  can  be  stabilized  globally  and  the  limit  cycle  will  be  eliminated. 
With  small  control  authority  the  region  of  attraction  of  the  origin  will  be  extended,  but  the  limit  cycle 
will  persist,  possibly  modified.  In  general,  one  cannot  claim  that  the  locally  stable  limit  cycle  will  be 
reduced  in  size,  even  though  this  it  typically  expected.  In  fact,  one  can  construct  examples  where  the 
control  action  that  extends  the  domain  of  attraction  of  the  origin  actually  increases  the  size  of  the 
locally  stable  limit  cycle. 

In  the  case  of  limited  control  authority,  when  the  elimination  of  the  limit  cycles  is  impossible,  a 
control  design  aimed  at  reducing  the  size  of  limit  cycles  usually  involves  approximate  methods  based 
on  harmonic  balance  (cf.  describing  function  methods)  where  the  controller  is  designed  using  the 
harmonic  balance  equation. 

In  paper  [55]  we  investigated  alternative  methods  of  reducing  the  size  of  large  limit  cycles  using 
control  in  the  case  of  limited  authority.  We  consider  a  planar  nonlinear  oscillatory  system  with 
bounded  control.  It  is  shown  that  under  weak  assumptions  one  can  identify  the  controller  that  shrinks 
the  existing  limit  cycles  and  maximizes  the  domain  of  attraction  of  the  origin.  Away  from  a  small 
neighborhood  of  the  origin,  the  optimal  control  action  results  in  switching  between  the  positive  and 
negative  saturation  values.  The  switching  happens  on  the  curve  where  transverse  linear  controllability 
is  lost  and  a  collection  of  vector  fields  drilling  the  local  control  properties  of  the  system  change 
orientation.We  briefly  compared  the  optimal  approach  to  controllers  suggested  by  a  number  of  well 
known  approximate  techniques. 

The  control  concepts  are  illustrated  using  the  Greitzer  compressor  surge  model.  In  the  situation 
when  the  throttle  is  set  so  that  the  uncontrolled  system  exhibits  a  surge  cycle,  a  control  via  small 
throttle  variation  can  reduce  or  eliminate  the  surge  cycle.  It  is  also  shown  that  in  the  case  of  the  surge 
model  the  controller  that  locally  minimizes  the  rate  of  divergence  from  the  origin  (in  the  given  metric) 
is  not  optimal  in  the  sense  of  convergence  to  a  smallest  possible  periodic  orbit. 

We  illustrated  the  results  on  Greitzer  compressor  surge  model.  Figure  3.1  shows  the  uncontrolled 
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Figure  3.1:  Uncontrolled  and  controlled  surge  cycle 


and  controlled  surge  limit  cycles  when  the  optimal  control  is  used. 


3.2  Adaptive  control  of  combustion  instability 

Adaptation  of  control  algorithm  parameters  to  optimize  the  performance  was  demostrated  at  UTRC 
in  combustion  instability  and  in  diffuser  flow  separation  control.  AFOSR-sponsored  research  played 
an  important  role  in  this  development.  Details  are  presented  in  papers  [cl0,c5,c7j4]. 

Experiments  and  model-based  analysis  of  combustion  instability  determined  that  pressure  mea¬ 
surement  and  a  simple  phase-shifting  algorithm  with  an  appropriately  chosen  control  phase  is  sufficient 
for  suppression  of  oscillations,  given  enough  control  authority  ([81,  82,  69,  94,  88,  96,  63,  70,  53,  54]). 

The  goal  of  using  active  combustion  instability  control  on  a  gas  turbine  engine  is  to  keep  pressure 
oscillations  at  an  acceptable  level  over  a  large  range  of  operating  conditions,  without  a  perfect  model 
of  the  process  dynamics.  Model  based  analysis  determined  that  minimum  information  needed  to  cal¬ 
culate  the  best  control  phase  for  the  combustion  instability  control  requires  estimation  of  parameters, 
including  transport  delay  (or  at  least  the  corresponding  phase  shift),  that  are  hard  to  obtain  from 
pressure  measurements  alone  [70,  53,  54].  Even  if  this  difficulty  could  be  circumvented,  the  sensitivity 
to  modeling  errors  was  likely  to  be  high. 

Need  for  developing  an  algorithm  that  would  allow  finding  the  best  phase  with  minimum  amount 
of  a  priori  information  that  would  work  over  large  range  of  operating  conditions  and  with  minimum 
model  assumptions  seemed  apparent.  The  operating  conditions  include  fast  engine  acceleration  and 
deceleration  transients. 

The  algorithms  developed  at  UTRC  and  presented  in  [56,  59,  60,  109]  allow  fast  automatic  tuning 
of  control  parameters  and  (under  reasonable  assumptions)  can  be  guaranteed  not  to  amplify  the 
oscillations.  Two  extremum-seeking  algorithms  [98,  99,  78]  have  been  selected  for  the  tuning  of  the 
control  phase.  Model  analysis  indicated  that  the  combustion  process  can  include  a  simple  pressure 
magnitude  dynamics  and  fast  stable  dynamics  describing  other  elements  of  combustion  process  and 
actuator  dynamics.  When  an  extremum-seeking  algorithm  is  applied  to  such  a  model,  it  is  possible 
to  guarantee  both  convergence  of  the  control  parameter  to  the  optimal  value  and  the  stability  of  the 
overall  system.  The  analysis  uses  only  the  fact  that  the  variable  to  be  minimized  (in  our  case  the 
pressure  magnitude)  has  an  almost  static  dependence  on  the  control  parameter.  The  stability  proof  is 
possible  if  one  exploits  the  time  scale  separation  between  the  slow  dynamics  of  the  control  parameter 
update  law  and  the  fast  dynamics  of  pressure  magnitude.  There  was  a  technical  challenge  in  finding  a 
phase  shifting  mechanism  that  works  for  pressure  oscillations  varying  over  large  frequency  range  and 
in  developing  a  reliable  pressure  magnitude  detection  mechanism.  This  challenge  has  been  overcome 
with  a  frequency  tracking  observer  algorithm  (based  on  Extended  Kalman  Filter  described  in  [83]) 


that  allows  for  fast  and  reliable  estimation  of  the  in-phase  and  quadrature  components  of  the  bulk 
pressure  mode  over  large  range  of  frequencies. 

Performance  specifications  for  extremum-seeking  algorithms  have  been  defined  for  algorithm  ini¬ 
tialization  transients  and  engine  acceleration  transients.  When  initialized  with  a  phase  corresponding 
to  amplification  of  oscillations,  the  algorithms  should  quickly  produce  and  maintain  phases  corre¬ 
sponding  to  suppression  of  the  oscillations.  In  the  engine  acceleration  transients  the  algorithms  should 
be  able  to  suppress  oscillations  relative  to  uncontrolled  levels. 

Experiments  with  an  extremum-seeking  adaptive  scheme  were  conducted  on  a  4MW  single  nozzle 
rig  at  United  Technologies  in  August  1998.  The  control  gain  was  fixed  and  only  the  control  phase  has 
been  updated.  One  of  the  algorithms  relied  on  estimation  of  the  derivative  of  the  pressure  magnitude 
with  respect  to  control  phase  by  introducing  a  small  sinusoidal  variation  in  the  control  phase  and 
measuring  the  response  of  the  pressure  magnitude  at  the  corresponding  frequency.  The  mean  control 
phase  was  incremented  by  an  amount  proportional  to  the  estimated  derivative. 

Initialization  transients  were  introduced  where  the  initial  control  phase  varied  significantly  from  the 
optimal  one.  The  algorithms  behaved  very  well  at  high  power  condition  (small  pressure  oscillations 
and  medium  level  of  broad-band  disturbance)  and  reasonably  well  at  low  power  conditions  (large 
pressure  oscillations  and  large  level  of  broad-band  disturbance).  Once  reaching  a  neighborhood  of 
the  optimal  value,  the  control  phase  usually  stayed  in  a  reasonably  small  neighborhood  of  that  value, 
rarely  produced  coP'rol  phases  corresponding  to  level  higher  than  uncontrolled  levels,  and  always 
provided  better  average  pressure  oscillations  levels  than  uncontrolled  levels.  Figure  3.2  shows  typical 
performance  of  the  algorithms.  Some  control  experiments  were  conducted  with  the  sole  purpose  to 
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Figure  3.2:  Extremum-seeking  experiments:  initialization  transients.  The  initial  control  phase  is 
strongly  destabilizing.  An  optimal  stabilizing  phase  is  reached  within  three  seconds. 

identify  the  pressure  magnitude  dynamics  and  the  noise  to  obtain  a  simple  model  of  the  slow  part  of 
the  combustion  process.  This  research  was  supported  by  the  AFOSR  grant.  The  model  obtained  at 
UTRC  by  Kartik  Ariyur  working  under  direction  of  Andrzej  Banaszuk  was  used  to  study  the  adaptive 
algorithms  off-line  [60]  and  simulate  fast  engine  transients  that  could  not  be  introduced  in  the  rig 
[56,  59]. 

A  detailed  model-based  analysis  of  the  adaptive  algorithm  based  on  averaging  was  conducted  by 
Miroslav  Krstic  from  UCSD  and  his  student,  Kartik  Ariyur.  This  analysis  provides  two  linear  models, 
one  for  tracking  reference  changes  and  the  other  for  sensitivity  to  noise,  which  offer  insight  into  how 
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different  parameters  influence  the  performance  [79,  60].  Details  are  presented  in  papers  [cl0,c5,c7j4]. 

3.3  Control  of  vortex  motion  and  chaotic  advection 

3.3.1  Introduction 

This  section  is  an  extended  abstract  of  the  work  described  in  detail  in  [c9J3]. 

In  the  last  two  decades,  much  progress  has  been  made  in  employing  dynamical  systems  theory 
to  laminar  mixing  [1,  2,  46,  48].  In  much  of  this  work  the  emphasis  has  been  on  understanding  the 
mechanisms  of  mixing  in  laminar  flows.  For  example  Melnikov  theory  and  associated  lobe  dynamics 
have  been  used  to  study  transport  and  mixing  in  and  out  of  the  recirculation  bubble  of  the  vortex 
Batchelor  couple  in  a  pioneering  study  of  Rom  Kedar  et  al  [47].  Those  mechanisms  are  now  well- 
understood  in  two-dimensional  flows  and  partially  in  three-dimensional  flows.  In  AFOSR-sponsored 
paper  [44]  we  pursued  a  related  question  of  optimizing  and  controlling  transport  and  mixing  in  two- 
dimensional  flows.  Akin  to  first  chaotic  advection  studies  we  chose  a  simple  vortex  model  -  a  single 
vortex  in  a  corner  subject  to  a  potential  field.  The  motivation  for  the  present  study  comes  from 
studying  flows  in  a  combustor.  (Vortex  models  have  been  used  before  for  this  purpose  [31],  but 
without  forcing.) 

In  the  absence  of  time-dependent  forcing  (introduced  by  modulating  the  potential  field)  the  vortex 
wass  either  at  a  stable  equilibrium  position  or  was  moving  on  a  periodic  trajectory  in  the  plane.  The 
first  question  that  we  asked  was:  if  we  allow  for  the  arbitrary  time-modulation  of  the  potential  field, 
can  we  move  the  vortex  from  an  arbitrary  initial  position  to  an  arbitrary  final  position  in  the  plane? 
The  answer  to  this  question  is  shown  to  be  positive  using  the  transformation  to  the  so-called  flat 
coordinates  [43,  18].  Having  this  result,  we  could  proceed  to  search  over  all  the  trajectories  of  the 
vortex  in  a  bounded  domain  of  the  plane  and  determine  the  optimal  one  with  respect  to  a  suitable 
measure.  We  chose  this  measure  to  be  the  flux  through  the  heteroclinic  orbit  connecting  two  stagnation 
points  of  the  partivcle  dynamics  (called  separatrix),  thus  linking  control  theory  and  chaotic  advection 
theory.  Once  the  optimal  trajectory  was  found,  we  stabilized  it  using  a  feedback  law.  The  search 
for  the  optimal  trajectory  was  pursued  here  using  the  numerical  simplex  method.  Finally  we  showed 
that  the  problem  presented  here  could  be  put  in  the  context  of  Pontryagin  maximum  principle  which 
provides  a  necessary  condition  for  an  optimal  trajectory  without  utilizing  the  transformation  to  flat 
coordinates. 

The  results  apply  to  the  case  of  single  vortex  motion  in  an  arbitrary  domain  driven  by  a  potential 
velocity  field  with  streamfunction  'J'c. 

For  more  details  see  [c9  j3]. 


3.3.2  Flattness-based  control  of  vortex  motion  in  a  combustor  recirculation  region 
The  equations  of  motion  for  the  vortex  are  given  by 


X*i  — 
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dyv  4 
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where  e(t)  is  the  control  input.  The  streamfunction  of  the  control  field  can  again  be  chosen  as  a  flat 
coordinate,  z\  =  ^c.  The  other  coordinate,  z 2  is  chosen  such  that  i\  —  z\ 2.  Thus, 


Z2  =  {^c«,}  = 


d9c  99 
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where  {tfc,  '&}  is  the  usual  notation  for  a  Poisson  bracket  of  functions  $c,  $.  It  is  easy  to  see  that 

i2  =  {{*c,  ¥},  *  }  +  «(*){{*« *  },  *  c}- 

Thus,  the  motion  of  a  single  vortex  in  a  container  of  arbitrary  shape  under  the  influence  of  a 
potential  field  with  a  streamfunction  is  locally  controllable  if 

{{*c,tf},tfc}^0  (3.1) 

on  the  domain.  (Another  generalization  holds  for  the  case  of  n  vortices  under  the  influence  of  n 
potential  fields.  See  [44]  for  details.  ) 

In  the  corner  problem  the  instantaneous  flow  is  determined  by  the  vortex  position.  Under  periodic 
conditions,  i.e.  at  natural  or  controlled  periodic  vortex  motion,  the  fluid  motion  may  be  characterized 
by  the  Poincare  map  x*  =  4>r(x),  where  x  represents  the  initial  position  of  the  fluid  particle  and  x*  the 
position  one  period  Tper  later.  At  small  amplitudes,  the  Poincare  map  has  (at  least)  three  fixed  points, 
one  at  the  origin,  one  xs  on  the  x-axes  and  one  xu  on  the  y-axis  (see  Fig.  3.3).  The  fixed  point  xs  {xu) 
has  a  stable  (unstable)  manifold  Wa(x4)  (W„(xu))  extending  into  the  domain.  The  manifolds  intersect 
each  other  infinitely  many  times  near  the  fixed  points  [23].  Fixed  points  and  invariant  manifolds  can 
be  shown  to  be  symmetric  with  respect  to  the  x  =  y  bisector,  if  the  vortex  motion  has  the  same 
symmetry  and  if  the  initial  position  of  the  vortex  is  on  the  x  =  y  line.  We  define  a  curve  C  which 
consists  of  the  sections  of  the  invariant  manifolds  from  the  fixed  points  to  the  first  intersection  point 
(Fig.  3.3).  For  the  considered  parameters,  the  curve  is  nearly  circular.  The  recirculation  region  R 
shall  be  defined  by  the  region  bounded  by  the  x-  and  y-axis  and  C.  In  the  limit  of  small  amplitudes 
R,  the  curve  C  converges  to  the  steady  state  separatrix,  except  in  an  arbitrary  small  neighborhood 
around  the  fixed  point.  In  the  sequel,  vortex  motion  is  assumed  to  be  periodic. 

Mixing  is  associated  with  the  behaviour  of  an  ensemble  of  fluid  particles  over  a  finite  period  of 
time.  A  necessary  condition  for  mixing  between  the  free  stream  and  the  recirculation  region  is  the 
fluid  exchange  across  C,  realizing  that  the  net  flux  vanishes  for  incompressible  flow.  The  instantaneous 
rate  of  fluid  exchange  is  quantified  by 

$inst  =  jds  |  Un| 

where  ds  represents  an  arc  element  of  C  and  un  the  normal  component  of  the  velocity.  A  suitable 
mixing  measure  $  may  be  defined  by  the  flux  averaged  over  one  period  and  normalized  with  the  area 

i*i. 

Tp*r 

* = ^  h  h  K|-  M 

Numerous  measures  for  mixing  can  be  considered  depending  on  the  application.  Examples  are  (1)  the 
flux  according  to  Eqn.  (3.2),  (2)  the  fluid  discharge  during  one  period,  (3)  the  size  of  the  mixing  region, 
(4)  the  average  residence  time  of  the  transient  fluid  particles  in  the  mixing  region,  (5)  the  average 
strain  rate,  since  the  strain  is  related  to  stretching  of  material  lines  and  surfaces,  (6)  the  growth  rate  of 
material  lines  and  surfaces,  (7)  the  average  divergence  rate  of  neighboring  fluid  particles  (Kolmogorov- 
Sinai  entropy),  (9)  any  measure  for  the  degree  of  finite-time  homogenization  of  an  ensemble  of  fluids 
particles.  In  [44]  we  restricted  to  a  few  of  these  quantities  which  can  be  expected  to  be  of  importance 
in  the  combustor.  In  particular  case  of  the  comer  flow  the  lobe  of  the  stable  manifold  Wt  contains 
the  fluid  which  leaves  the  recirculation  region  R  during  the  period  described  by  the  Poincar£  map 
(Figure  3.3).  This  fluid  is  bounded  by  the  first  intersection  point,  the  W4-lobe  and  the  initial  Wu 
section.  The  lobe  of  the  unstable  manifold  Wu  contains  the  fluid  which  enters  R  during  the  same 
interval.  This  fluid  is  bounded  by  the  first  intersection  point,  the  Wu-lobe  and  the  initial  Ws  section. 
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Hence,  the  intersection  of  both  areas  identifies  the  fluid  which  enters  and  leaves  during  the  considered 
time  interval,  i.e.  the  area  between  invariant  manifolds  connecting  intersection  points  2  and  3.  Thus, 
the  total  fluid  exchange  Tper  |-R|$  during  one  period  is  represented  by  the  areas  bounded  by  C  and  both 
lobes,  W3,WU,  counting  the  overlap  area  twice.  The  region  around  the  vortex  which  is  embraced  by 
both  lobes  neither  leaves  nor  enters  the  recirculation  zones  during  the  considered  period.  The  embraced 
region  can  roughly  be  identified  as  the  trapped  region  neglecting  convection  effects  due  to  the  thin 
higher-order  lobes  near  the  axes  (not  resolved  by  present  computations).  Thus,  the  Poincare  map 
allows  to  draw  many  quantitative  and  qualitative  conclusions  about  the  mixing  measures  discussed  in 

[44]- 


Figure  3.3:  Invariant  manifolds  for  the  Poincare  map  of  the  controlled  particles  motion 

Figure  3.4  illustrates  the  concept  of  using  vortex  model  and  methods  from  control  theory  and 
dynamical  systems  theory  for  control  of  mixing  in  the  corner  flow  example. 

For  more  details  see  [c9J3]. 

3.3.3  Observer  based  control  of  vortex  motion  in  a  combustor  recirculation  region 

The  tracking  control  developed  above  postulates  the  availability  of  vortex  coordinates  for  feedback. 
To  become  part  of  a  plausible  implementation  framework  it  therefore  needs  to  be  complemented  by  a 
dynamic  observer  that  translates  some  feasible  sensor  readings  to  state  estimates.  This  section  briefly 
summarizes  aspects  of  an  observer  based  design.  More  details  are  provided  in  [ell  j5]. 

The  assumed  measurements  are  of  the  s-direction  fluid  velocity,  denoted  u,  at  the  normalized  wall 
point1  (1,0).  The  following  is  an  explicit  expression  for  u 

U  =  dy(*av  +  (1  +  e)^c)|x=l,y=0  =  0.5(1  +  e)  —  ^2  +  y2  +  1)2  _  4^2 

The  component  0.5(1  -I-  e)a  adds  no  information  and  will  be  ignored. 

1  Since  oar  model  does  not  account  for  viscosity,  the  idealization  of  velocity  measurement  at  the  wadi  is  admissible. 
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Analysis  of  mixing  using  control  and  dynamical  systems  theory  tools 


Example:  control  of  mixing  In  a  combustor 
Control  theory  applied  to  vortex  dynamics  model : 

differential  flatness  of  the  rortex  dynamics  allows  to  create 
and  stabilize  periodic  vortex  motion 

Dynamical  systems  theory  applied  to  particle  dynamics  model: 

geometric  properties  of  Invariant  manifold  of  Poincare  map 
determine  particle  mixing  properties 


dxjdt  *  ffx+t  )  vortex  dynamics 

dx/dt  -//x„  x,e  )  particle  dynamics 
xr-  vector  of  vortex  position 
xf-  vector  of  particle  position 
e  -  control 


Figure  3.4:  Illustration  of  control  of  mixing  in  the  corner  flow. 


Straightforward  analysis  (left  out)  verifies  that  local  linearizations  of  our  system  are  observable 
almost  everywhere,  relative  to  the  output  u,  and  that  the  set  (curve)  formed  by  points  of  lost  observ¬ 
ability  are  not  invariant.  Nonetheless,  we  chose  not  to  attempt  to  design  a  global  observer.  This  was 
motivated  by  the  level  of  nonlinearity  of  our  system  (including  the  vortex  model  and  (3.3)  and  the  fact 
that  a  global  observer  would  require  sensitive,  state-estimate  dependent  gain  scheduling.  For  exam¬ 
ple,  an  observer  based  on  the  canonical  local  coordinates  (u,  u)  [76,  58],  would  require  back-and- forth 
coordinate  change  ( xv,yv )  <4  (u,  u).  Instead  we  opted  for  a  local  but  more  robust,  fixed  structure 
observer  that  will  cover  a  neighborhood  about  the  fixed  point  to  which  designated  tracking  orbits  will 
be  restricted.  This  observer  will  be  used  in  conjunction  with  the  state  feedback  policy  proposed  in 
previous  section.  It  will  be  complemented  by  a  dissipative  feedback  control,  used  to  drive  the  vortex 
into  the  designated  neighborhood.  The  latter  does  not  require  a  detailed  state  estimate,  and  makes 
do  with  direct  feedback  of  the  measurement  signal. 


3.3.4  A  Dissipative  Constant  Structure  Observer  for  Small  Deviations 

To  match  the  tracking  feedback,  the  suggested  observer  is  written  in  the  flat  coordinate  system,  and 
is  thus  form 

Z\  =  i2  +  Ci(fi  — u)  .  #  ■  (3_4) 

'z2  =  p  +  eq  +  C2(u  -  «) 


where  a  hat  indicates  an  estimated  variable  and  where  (i  are  correction  terms.  The  simplest  option  is 
that  of  using  a  constant  linear  gain  - 
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where  <5  indicates  an  estimation  error.  Using  linear,  small  deviations  approximations,  the  dynamics  of 
state  estimation  errors  are  then  governed  by  - 


6z\ 

-Lidziu 

1  —  L\dZ2xi 

f  _  \ 

Sz\ 

dt 

Sz2 

—  jr(l  —  z|)  —  LzdZlu 

L 

2z2  ^1  +  e  —  — J  —  L2dZ2u 

Sz2 

and  the  selection  of  gains  Li  should  stabilize  this  system. 

An  intuitive  guidance  for  such  selection  is  based  on  assuring  stability  of  the  frozen  time  systems: 
it  is  an  obvious  fact  that  systems  of  the  form 
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Sz2 

with  (constant)  b,  c  >  0,  a  >  0  and  d  >  0,  generates  stable  dynamics.  Direct  evaluation  reveals  that 
the  signs  of  both  partial  derivatives  dZiu  are  both  positive  in  a  neighborhood  of  the  rest  point  (1,1). 
Following  the  frozen  time  intuition,  this  suggests  using  L\  =  0  and  Lz  >  0. 

Since  stability  of  frozen  time  system  neither  implies,  nor  is  it  implied  by  stability  of  the  actual, 
time  varying  system,  the  latter  has  to  be  analyzed.  To  the  very  least,  such  an  analysis  needs  to  be 
applied  to  a  periodic  linearization,  in  a  neighborhood  of  the  designated  reference  trajectory.  A  useful 
tool  here  is  to  replace  that  periodic  system  by  a  linear  time  invariant  model  for  approximate  dynamic 
phasors  of  the  error  signal .  The  term  dynamic  phasors  is  used  in  reference  to  the  time  varying  Fourier 
coefficients  of  a  trajectory  that  is  dominated  by  a  known  period  T  - 


<  x  >k  (t)  =  ±  J°  x(t  +  r)e~^t+^dr  (3.8) 

Differential  equations  governing  the  evolution  of  dynamic  phasors  can  be  formally  derived  from  the 
differential  equation  for  the  original  time  trajectory,  using  the  equality  - 


^<x>k{t)  =  -ju  <x>k  (<)+  <  x  >k  ( t )  (3.9) 

where  to  =  2 tt/T.  Approximate  models  are  based  on  truncated  Fourier  series  expansions.  Dynamic 
phasors  of  a  T-periodic  trajectory  are  constant,  and  a  desired  periodic  trajectory  is  thus  represented 
by  a  desired  set  point  in  phasor  representation.  The  linearized  dynamic  phasor  model,  about  a  desired 
set  point  is  thus  time  invariant.  Local  stability  analysis  and  guidance  for  parameter  selections  thus 
reduces  to  eigenvalue  analysis  in  the  respective  model.  In  the  present  case,  such  analysis  (left  out 
here)  reveals  that  stability  is  attainable  for  reference  trajectories  with  a  period  of  u  =  1,  dominated 
by  their  first  harmonic  (in  the  flat  coordinates),  centered  at  the  rest  point  and  with  a  radius  of  ~  0.5. 


3.3.5  A  Simple  Dissipative  Measurement  Feedback 

The  dynamic  observer  suggested  above  is  restricted  to  periodic  orbits  in  a  neighborhood  of  the  rest 
point.  An  additional  control  mechanism  is  therefore  needed  to  drive  the  vortex  into  such  a  neighbor¬ 
hood.  In  the  absence  of  an  observer,  such  a  mechanism  cannot  require  detailed  state  information.  The 
tool  of  choice  here  is  the  equation  - 


+  (1  +  e)*c)  =  (3-10) 

By  (3.10),  any  control  selection  that  forces  e'I'c  <  0,  is  dissipative,  and  will  drive  the  vortex  towards 
the  rest  point.  Moreover,  since  'I'c  is  close  to  periodic  under  low  gain  actuation,  the  same  holds  if  i 
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is  selected  to  be  (close  to)  zero  mean,  with  e($c  -  tfc)  <  0,  where  is  the  average  over  a  period. 
The  advantage  of  a  zero-mean  e  is  that  it  prevents  a  drift  in  e,  and  results  with  a  practically  feasible 
bounded  actuation  signal. 

Figure  3.5  shows  plots  of  \I,C  -^>c  and  of  a  low-pass  version  of  u  -u;  both  trajectories  are  associated 
with  a  vortex  trajectory  that  slowly  spirals  towards  the  rest  point.  The  phase  match  is  obvious  and 
a  path  to  a  stabilizing  control  is  set.  Figure  3.6  depicts  the  stabilizing  effect  of  a  low  gain  negative 
integral  feedback  e  =  — 6(u  —  u)  —  pe,  where  S  >  0  is  a  small  gain  and  p  >  0  is  a  small  forgetting 
factor”  that  is  used  to  eliminate  long  term  effects  of  deviations  from  the  ideal  zero-mean  condition. 


Figure  3.5:  Plots  of  zx  -  z\  (solid)  and  a  (low  pass  filtered)  u-u  (dashed)  under  low  gain  actuation. 


Figure  3.6:  Phase  portrait  of  vortex  motion  under  the  simple  dissipative  feedback 

One  practical  aspect  that  needs  to  be  resolved  in  implementing  the  dissipative  control  suggested 
above  is  a  means  of  estimating  the  period  of  motion  during  weakly  actuated,  “large  deviations” .  Such 
estimates  are  needed  in  order  to  extract  tZ  from  velocity  measurements.  Additionally,  an  estimate  of 
the  total  energy  stored  (based  on  the  combined  Hamiltonian  'Pmt)  -I-  e^c)  is  needed,  in  order  to  be 
able  to  determine  whether  the  vortex  entered  (or  left)  the  domain  of  the  dynamic  observer.  Direct 
inspection  reveals  that  (in  the  weakly  actuated  system),  there  is  a  monotonous,  easily  tabulated 
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relationship  between  both  these  quantities  and  the  minimum  over  a  period  of  the  measured  velocity  u 
as  depicted  in  Figure  3.7.  Since  period  variations  (in  the  normalized  coordinates)  are  roughly  bounded 
between  [27r,37r],  the  estimation  of  both  quantities  -  without  a  priori  information  -  requires  less  then 
two  periods. 


Figure  3.7:  Plots  of  the  approximate  value  of  z\  -  when  crossing  the  level  z<i  —  0  (left)  -  and  of  the 
approximate  period  (right)  in  unactuated  dynamics,  as  functions  of  the  minimal  absolute  value  over 
a  period  of  the  measured  fluid  velocity 


3.3.6  A  Global  Observer  Based  Tracking  Controller 

A  global  controller  is  formed  by  patching  the  large  deviations  dissipative  controller,  as  described 
right  above,  with  a  controller  that  combines  the  state-feedback  tracking  controller  and  the  dynamic 
state  observer,  in  a  neighborhood  of  the  prescribed  vortex  trajectory.  As  is  usually  the  case,  safe 
switching  between  different  controllers  is  a  delicate  job.  Critical  issues  include  stability  under  the 
switching  procedure  and  using  hysteretic  bands  in  the  switching  region,  to  avoid  chatter.  Details  of 
these  mechanisms  me  not  specific  to  the  focal  points  of  this  article  and  are  left  out.  A  demonstration 
of  a  response  of  the  bi-modal  closed  loop  system  is  depicted  in  Figure  3.8,  where  the  controller  -  in 
dissipative  mode  -  drives  the  vortex  from  the  far  field  to  the  desired  neighborhood,  and  then  switches 
to  steady  state  tracking  of  the  designated  trajectory. 

3.3.7  Remarks  On  Generic  Aspects 

As  some  details  of  the  developments  described  in  this  section  may  appear  to  be  peculiar  to  the  specific 
system  at  hand,  it  is  worth  highlighting  those  aspects  that  are  likely  to  be  generic  to  similar  control 
problems  involving  low  order  vortex  models. 

•  The  first  comment  concerns  the  emphasis  on  a  simple  structure  /  constant  gain  observer.  The 
systems  under  consideration  are  characterized  by  the  facts  that  detailed  models  are  both  distributed 
and  highly  nonlinear,  and  that  nonlinearity  is  a  dominating  aspect  of  even  the  simplest,  lumped,  low 
order  models,  used  in  control  design.  Absent  reliable  modeling  error  bounds,  avoiding  the  use  of  a 
continuum  (or  a  very  large  number)  of  local  linearizations  and  elaborate,  state  estimate  dependent 
coordinate  change  mappings,  seems  to  be  a  matter  of  common  sense  prudence.  Our  example  illustrates 
the  fact  that  the  domain  of  validity  of  a  fixed  structure  observer  may  be  large  enough  to  cover  the 
desired  region  of  operation.  In  other  systems,  switching  between  a  small  number  of  fixed  observers  may 
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Figure  3.8:  Closed  loop  tracking  of  a  vortex  orbit  reference,  asterisks:  tracking  reference  orbit;  bold: 
actual  vortex  trajectory 

be  necessary  to  cover  the  area  of  interest,  but.  the  objective  of  avoiding  a  continuous  gain  scheduling 
may  still  be  attainable. 

•  When  both  free  and  controlled  dynamics  are  periodic,  direct  readings  of  “over  a  period”  properties  of 
the  sensed  signal  may  provide  valuable  information.  In  particular,  it  may  provide  the  triggers  needed 
to  switch  between  different  modes  of  a  multi  modal  controller,  or  even  for  direct  feedback  stabilization, 
as  demonstrated  by  the  present  example. 

•  Continuing  the  focus  on  periodicity,  linearized  dynamic  phasor  models  provide  a  useful  tool  for 
stability  analysis  and  design  parameter  assignments. 

•  Its  simplicity  notwithstanding,  a  very  interesting  observation  concerns  the  dissipative  mode  of  the 
suggested  controller.  While  standard  dissipative  controllers  are  based  on  adding  or  increasing  some 
natural  form  of  damping  (e.g.,  current  resistance  or  mechanical  friction),  our  vortex  model  is  inherently 
lossless  in  steady  operation  and  control  is  not  used  to  add  damping  in  the  above  sense!  Here  we  use 
an  altogether  different  mechanism,  based  on  the  fact  that  the  total  stored  energy  is  the  sum  of  two 
distinct  Hamiltonians  {'t’mv  and  $c)  and  that  control  actuation  is  used  to  vary  the  relative  weight 
of  each  of  the  two  Hamiltonians  in  the  combined  energy  term.  In  words,  dissipation  is  achieved  by 
varying  the  relative  weight  of  in  opposite  phase  to  the  variations  in  'Fc. 


3.4  Diffuser  pressure  recovery  control 

UTRC  has  range  of  flow  control  programs,  including  DARPA  Micro  Adaptive  Flow  Control  program 
aimed  at  helicopter  retreating  blade  stall  control.  Recently,  dynamical  systems  and  control  theory 
methods  were  applied  to  control  of  flow  separation  in  diffusers  within  the  Dynamic  Separation  Control 
project.  This  effort  was  partially  sponsored  by  AFOSR.  Details  are  decribed  in  paper  [c6]. 

Aggressively  open  diffusers  are  subject  to  flow  separation  and  hence  loss  of  performance  measured 
by  the  pressure  recovery  coefficient.  This  coefficient  is  defined  by  Cp  =  t^j,  where  A p  =  pi  —  pi 
represents  the  static  pressure  difference  between  inlet  and  exit,  p  the  density  and  U\  the  mass-  and  time- 
averaged  inlet  velocity.  The  subsonic  diffuser  performance  measured  by  pressure  recovery  coefficient 
can  be  improved  by  increasing  mixing  of  the  high  momentum  free  stream  and  low  momentum  fluid 
at  the  wall.  Separation  control  in  a  diffuser  at  large  angles  of  expansion  was  acchieved  at  UTRC 
experiments  described  in  [38, 42, 12]  via  periodic  forcing  using  a  synthetic  jet  actuator  located  usptream 


of  the  separation  point.  (Previously  such  concepts  were  used  to  control  of  airfoil  separation  [49]). 
Figure  3.9  shows  the  experimental  setup.  In  particular,  in  [12]  UTRC  team  and  Brianno  Coller 


reference  station 
(monitor  inlet  velocity  and  static  pressure) 


q  dynamic  pressure 

P(0  Ap(t) 
o  #  o  #o  o 


forcing 
location 


Figure  3.9:  Flow  facility  schematic 

has  studied  diffuser  pressure  recovery  control  using  models  and  experiments.  The  new  concepts  in 
this  study  was  to  use  control  to  excite  beneficial  non-equilibrium  dynamics  in  the  separating  shear 
layer.  Acoustic  forcing  at  the  beginning  of  the  expanding  section  of  the  diffuser  resulted  significant 
performance  enhancements  for  certain  actuation  frequencies  as  shown  in  Figure  3.10.  Brianno  Coller 


Figure  3.10:  Pressure  recovery  coefficient  in  terms  of  the  frequency  at  the  diffuser  exit  plane  for  the 
forced  separated  Bow,  showing  sensitivity  of  performance  to  frequency. 

(in  collaboration  with  UTRC  reseachers)  has  developed  a  Lagrangian  point  vortex  simulation  [33,  35] 
to  predict  pressure  recovery  enhancement  in  a  planar  diffuser  as  a  result  of  actuation  [12].  Numerical 
results  compared  well  with  experiment.  Qualitative  evolution  and  interaction  of  the  controlled  vortices 
near  the  wall  are  viana.li7.ftd  using  smoke  and  are  also  found  to  match  the  model  simulation  results 
well. 
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Figure  3.11:  Vortex  formation  of  a  diffuser  Bow  forced  with  a  Strouhal  frequency  of  St  =  0.3. 


For  more  information,  please  refer  to  paper  [c6]. 

3.5  Control  of  mixing  in  jets  in  cross-flow 

3.5.1  Introduction 

Mixing  in  various  places  in  combustor  is  the  key  to  aeroengine  performance  and  durability  of  its 
parts.  The  temperature  pattern  factor  is  a  measure  of  nonuniformity  of  the  temperature  profile  of 
the  air  approaching  the  aeroengine  turbine.  Nonuniformity  creates  excess  oxidation  and  large  thermal 
stresses  in  the  turbine  guide  vanes  that  lead  to  destruction  of  the  meted  and  hence  to  reduced  turbine 
durability.  On  the  other  hand,  because  of  high  engine  thrust  requirements  in  military  aeroengines,  one 
needs  to  operate  at  high  turbine  inlet  mean  temperature  conditions.  This  leads  to  high  maintainance 
costs,  as  the  vanes  must  be  replaced  often.  The  temperature  profile  is  controlled  by  the  mixing  rate 
of  the  air  jets  injected  into  the  combustion  chamber  along  the  combustor.  Improved  mixing  of  the  air 
jets  with  the  cross-flow  carrying  hot  combustion  products  would  allow  a  more  uniform  temperature 
profile  of  the  air  at  the  turbine  exit,  with  reduced  peak  temperature  and  gradients,  while  maintaining 
the  same  mean  temperature,  and  hence  thrust.  In  this  way  an  increase  of  mixing  rates  of  jets  in 
cross-flow  using  control  would  improve  aeroengine  affordability. 

3.5.2  Description  of  Dynamic  Combustion  Enhancement  project 

The  AFOSR  sponsored  research  at  UTRC  in  2000  was  focused  on  control-oriented  modeling  and 
model-based  control  of  jets  in  cross  flow.  This  work  was  an  integral  part  of  larger  effort  in  Dynamic 
Combustion  Enhancement  (DCE),  mostly  sponsored  by  corporate  sources.  The  goal  of  DCE  project  in 
2000  was  to  demonstrate  in  experiment  a  significant  mixing  enhancement  in  jet  in  cross-flow  using  flow 
control,  or  prove  using  a  model  that  this  is  impossible  using  low  actuation  level.  The  AFOSR  project 
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Figure  3.12: 

was  focused  on  development  of  tools  for  determining  actuator  placement  and  control  laws  based  on 
reduced-order  models.  DCE  project  used  the  tools  designed  within  AFOSR  project  and  demonstrated 
mixing  enhancement  with  control  in  experiments. 

3.5.3  Physical  phenomena  observed  in  jets  in  cross-flow 

It  was  determined  from  the  literature  studies  done  in  January  and  February  2000  and  from  results 
of  UTRC  experiments  and  Large  Eddy  Simulations  (LES)  that  the  dynamics  of  the  jet  in  cross-flow 
is  dominated  by  two  types  of  coherent  structures  (see  Figure  3.13).  The  shear  layer  Kelvin-Helmoltz 
type  structures  occur  where  the  jet  bends  in  the  direction  of  the  cross-flow.  Downstream  of  the  bend 
the  dynamics  is  dominated  by  streamwise  counter-rotating  vortices.  This  counter-rotating  vortex  pair 
(CVP)  seems  to  be  the  main  mixing  mechanism  of  the  uncontrolled  jet.  It  is  aggreed  in  the  literature 
that  CVP  is  the  major  mixing  structure  in  the  jets  in  cross-flow.  The  dominance  of  the  coherent 
structures  is  a  strong  indication  that  reduced-order  modeling  is  a  viable  modeling  path.  We  would 
like  to  mention  that  even  though  the  jets  in  cross-flow  have  been  subject  of  extensive  studies,  the 
mechanism  of  creation  of  CVP  is  still  a  matter  of  controversy  and  active  research. 

Any  study  of  miring  requires  investigation  of  particle  motion.  In  experiments  this  is  done  using 
Planar  Laser  Induced  Fluorescence  (PLIF)  or  MIE-scattering.  For  instance,  in  PLIF  technique  fluo¬ 
rescent  particles  are  introduced  to  the  jet  flow  and  fast  speed  camera  takes  pictures  of  a  planar  sheet 
where  a  pulsed  laser  excites  the  fluorescence  of  the  acetone.  Figure  3.14  shows  a  typical  experimental 
setting.  In  this  way  instantaneous  measurements  of  the  concentration  of  passive  scalar  introduced 
in  the  jet  can  be  obtained  with  spatial  resolution  of  30  microns  and  temporal  resolution  of  30Hz. 
Averages  of  many  snapshots  allow  to  obtain  spatial  probability  distribution  of  the  passive  scalar  in 
the  measurement  plane.  The  probability  distribution  allows  to  calculate  one  of  the  commonly  accept- 
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Jet  in  cross-flow:  unsteady  fluid  structures  determine  mixing 

Counter-rotatingVortex  Pair 


Figure  3.13:  Coherent  structures  dominating  in  jets  in  cross-flow 


Experimental  measurement  of  passive  scalar  concentration 


Figure  3.14:  Measurement  of  passive  scalar  concentration  in  PLIF  experiment 


ed  measures  of  mixing  called  unmixedness,  defined  as  the  ratio  of  the  standard  deviation  of  passive 
acalar  concentration  to  its  mean  value.  Decreasing  of  the  unmixedness  on  a  plane  transverse  to  the 
jet  located  downstream  of  the  jet  bend  by  a  significant  factor  is  the  goal  of  DCE  project  at  UTRC. 

Figure  3.15  shows  results  of  passive  scalar  measurements  in  UTRC  experiments.  The  analysis 
of  passive  scalar  pictures  shows  the  effect  of  the  counter-rotating  vortex  pair  on  the  passive  scalar 
distribution. 
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Figure  3.15:  Experimentally  measured  passive  scalar  concentration  in  jet  in  cross  flow  experiments  at 
UTRC 

3.5.4  Mathematical  description  of  the  flow 

We  are  interested  in  modeling  a  strong  jet  entering  from  a  pipe  into  a  crossflow.  A  schematic  repre¬ 
sentation  of  such  a  jet  is  shown  in  Figure  3.13.  The  quantities  Uqc  and  f7jet  are  the  velocities  of  the 
undisturbed  crossflow  and  of  the  jet,  and  the  radius  of  the  jet  is  Rjet-  We  start  with  the  assumption 
that  for  a  strong  jet  in  crossflow,  the  vorticity  is  concentrated  on  the  surface  of  the  jet  [11].  Then  the 
jet  can  be  modeled  by  a  cylindrical  vortex  sheet  governed  by  the  three-dimensional  inviscid  vorticity 

equation,  i.e. 

—  +  vTVw  —  wrVv  =  0  (3-11) 

dt 
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where  v  and  u  are  the  velocity  and  vorticity  vectors,  and  w  =  Vxv.  In  addition,  the  velocity  vector 
v  satisfies  the  incompressible  continuity  equation  [6] 

VTv  =  0  (3-12) 


Using  Eq.  (3.12),  Eq.  (3.11)  can  be  further  written  in  conservation  form 

^  +  (V7*Tf  =  0  (3.13) 

at 

n  here  F  is  the  3x3  flux  matrix,  i.e. 

F  =  cjvt  —  (3-14) 


The  closure  relation  between  the  vorticity  and  velocity  is  given  by  the  vorticity-induction  equation 
(Biot-Savart  law)  [86].  The  velocity  v  can  be  expressed  as  the  sum  of  an  irrotational  component,  the 
incoming  crossflow  velocity  Uqo,  and  a  solenoidal  component,  the  vorticity-induction  component  vw 


v  =  Uoo  +  vw  ,  Vu 


s/v"* 


(3.15) 


Passive  scalar  transport  is  described  by  the  ad vect ion- diffusion  equation 


dCfe--  +  v(x,  t)rVc(x,  t )  -  /iAc(x,  t)  =  0. 


(3.16) 


For  the  flows  of  interest  the  advection  term  v(x,  t)TVc{x,t)  dominates  the  diffusion  term  /iAc(i,  t). 

The  passive  scalar  distribution  from  any  model  can  be  compared  with  experimental  data  obtained 
using  the  setup  shown  on  Figure  3.14. 

The  control  of  the  jet  can  be  achieved  by  manipulation  of  boundary  conditions.  For  instance, 
the  profile  of  the  flow  entering  the  cross-flow  from  the  pipe  can  be  assumed  to  be  a  given  function 
of  time.  The  mixing  performance  can  be  measured  by  various  functionals  depending  on  the  passive 
scalar  concentration  c(x,t).  For  instance,  the  instantaneous  unmixedness  (the  ratio  of  the  standard 
deviation  to  the  mean  value)  on  the  performance  measurement  plane  V  downstream  of  the  jet  bend. 

The  vorticity  formulation  of  Navier-Stokes  equation  is  an  infinite-dimensional  evolution  equation 
w’th  the  right-hand  side  containing  a  nonlinear  integral  operator.  The  full  infinite-dimensional  model 
of  jet  in  cross-flow  is  clearly  intractable  for  analysis  and  control  design.  Reduced  order  modeling  is 
a  necessity.  The  current  effort  in  control-orientted  reduced-order  modeling  of  jets  in  cross  flow  is 
described  in  the  following  sections. 


3.5.5  Requirements  for  the  control-oriented  models 

Since  the  goal  of  the  project  is  to  make  an  impact  on  Dynamic  Combustion  Enhancement  future  using 
methods  of  dynamical  systems  and  control  theory,  we  require  the  models  to  be  physically  relevant , 
suitable  for  control  design,  and  refinable.  Here  is  a  more  precise  description  of  the  metrics  for  the 
models. 


1.  Physically  relevant.  The  models  should  adequately  capture  evolution  and  interactions  of  the 
most  dominant  coherent  structures  observed  in  the  flow.  In  particular,  this  excludes  input-output 
black-box  models  that  do  not  use  explicit  physical  parameters. 

2.  Suitable  for  control  design.  Models  need  to  be  suitable  for  control  design  using  existing 
methods  of  control  theory,  not  just  control  simulations.  In  particular,  this  implies  strong  preference 
for  the  following  requirements. 
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(a)  Explicit  apperance  of  control  input  as  a  function  of  time,  not  just  availability  of  few  tunable  pa¬ 
rameters. 

(b)  Ability  to  explicitly  calculate  an  instantaneous  measure  of  mixing  performance. 

(c)  Continuous  finite  dimensioned  description. 

3.  Refinable.  The  models  should  have  a  structure  that  allows  one  to  obtain  a  more  accurate  s- 
patial  and  temporal  description  of  the  flow  when  more  computational  elements  are  used.  In  the  limit, 
the  models  should  approach  CFD  accuracy  of  description. 


The  Lagrangian  vortex  models  have  many  advantages  over  Eulerian  CFD  for  simulations  of  shear 
flows  at  high  Reynolds  numbers  [33,  33,  35].  The  low  order  vortex  models  used  effectively  for  dynam¬ 
ical  systems  studies  of  mixing.  However,  purely  Lagrangian  models  are  not  best  suited  for  traditional 
control  design  for  the  flows  dominated  by  convection.  The  main  problem  is  that  the  computational 
vortex  elements  have  to  be  created  and  removed  during  discrete  events  determined  by  the  computa¬ 
tional  algorithm.  This  makes  Lagrangian  models  hybrid.  Another  disadvantage  of  Lagrangian  vortex 
models  is  that  the  control  input  and  surfaces  on  which  the  performance  measure  is  evaluated  are 
usually  fixed  in  the  Eulerian  frame.  Thus,  in  most  cases  of  interest  the  Lagrangian  models  do  not 
satisfy  requirement  of  suitability  for  control  design.  Therefore,  we  deeded  to  pursue  an  alternative 
modeling  path.  The  concept  of  mixed  Eulerian-Lagrangian  models  wr  i  introduced  by  Bernd  Noack 
and  Andrzej  Banaszuk.  A  generalization  of  the  concept  for  the  three-dimensional  vortex  sheet  was 
proposed  Razvan  Florea  and  is  currently  under  development. 

The  model  validation  metrics  has  been  defined  as  well.  It  has  been  determined  that  the  models 
should  agree  with  experimental  data  with  respect  to  the  following  four  metrics: 

1.  Jet  centerline  trajectory. 

2.  Counter-rotating  vortex  pair  peak  vorticity  trajectories. 

3.  Circulation  assigned  to  peak  vorticity  trajectory. 

4.  Passive  scalar  distribution. 

All  of  these  physical  quantities  can  be  obtained  from  experiments.  The  easiest  ones  to  obtain  are  the 
jet  centerline  trajectory  and  concentration  of  passive  scalar.  The  other  two  quantities  require  lengthy 
measurements  of  velocity  fields  by  transversing  a  hotwire  probe  or  a  costly  Particle  Image  Velocimetry 
(PIV)  measurement. 

3.5.6  Hierarchy  of  control-oriented  models  for  jets  in  cross-flow 

Hierarchy  of  models  for  jets  in  cross-flow  created  at  UTRC  includes  the  following  models: 

1.  High  complexity  model:  3D  spectral  Direct  Numerical  Simulation  (DNS)  for  jets  in  cross-flow 
including  passive  scalar  transport  equation  to  trace  concentration  of  jet  particles  was  developed  and 
transitioned  to  UTRC  by  Thomas  Bewley  and  Peter  Blossey  from  UCSD.  Code  allows  overnight 
simulations  for  both  uncotrolled  and  controlled  cases  on  UNIX  workstations  and  fast  PCs.  DNS  is 
considered  physically  relevant  and  refinable,  but  not  directly  suitable  for  control  design  (millions  of 
states).  The  code  was  used  for  extraction  of  parameters  of  grey-box  type  models  and  for  validation  of 
low  complexity  models  and  control  concepts  derived  using  these  models.  After  a  suitable  performance 
measure  (penalizing  nonuniform  distribution  of  passive  scalar)  has  been  identified  from  studies  of  low 
complexity  models  we  will  attempt  optimization  of  control  laws  and  actuator  locations. 

2.  Medium  complexity  models:  3D  Eulerian/Lagrangian  vortex  sheet  model  has  been  introduced 
by  Razvan  Florea,  Bernd  Noack,  and  Andrzej  Banaszuk  and  is  currently  being  developed.  The  jet 
flow  is  modeled  by  a  cylindrical  vortex  sheet.  Starting  with  the  incompressible  vorticity  equation,  an 
integral  formulation  for  the  strength  of  the  vortex  sheet  in  a  moving  coordinate  system  is  derived. 
The  velocity  field  induced  by  the  vortex  sheet  is  given  by  the  Biot-Savart  law.  The  computational 
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mesh  parametrizing  the  jet  surface  can  move  along  planes  locally  transverse  to  the  jet  axis.  This 
characterizes  a  “Lagrangian”  part  of  the  model.  The  transverse  planes  are  either  fixed,  or  can  rotate 
with  the  jet  surface  but  cannot  move  along  the  jet  axis.  This  characterizes  an  Eulerian  part  of  the 
model.  The  model  is  still  under  development.  Boundary  condition  need  to  be  fixed  and  the  simulations 
results  need  to  be  validated.  If  the  model  is  passes  a  validation  test,  it  offers  a  great  potential  for 
control  design.  It  is  physically  relevant,  and  refinable.  The  model  uses  a  fixed  number  of  states  (of 
order  of  thousands).  With  the  run  time  in  minutes  it  is  more  suitable  for  control  design  than  DNS. 
Unfortunately,  control  techniques  applicable  for  models  with  thousands  of  states  are  still  limited. 

3.  Low  complexity  model:  parabolized  two-dimensional  two-vortex  model  has  been  developed  by 
Thomas  John,  Satish  Narayanan,  Andrzej  Banaszuk,  and  Shubhro  Ghosh. 

3.5.7  Control  of  jet  in  cross  flow  using  DNS  model 

We  have  performed  direct  simulations  of  the  jet  in  crossflow  to  evaluate  the  effect  of  forcing  on  mixing 
by  the  jet.  Both  forced  and  unforced  cases  were  simulated  at  a  Reynolds  number  Rej  =  Ujdj/u  = 
3000  and  a  jet  velocity  ratio  R  =  Uj/Uoo  =  6.  A  number  of  forced  cases  were  explored  in  shorter 
simulations,  and  two  forced  cases  were  chosen  for  long  simulations  to  provide  well-converged  statistics 
for  comparison  with  the  unforced  case.  Both  the  time-averaged  properties  of  the  jets  and  the  level 
of  unsteadiness  in  mixing  are  evaluated  using  a  number  of  metrics,  including  jet  trajectory,  spreading 
rate,  entrainment  and  mixedness.  The  forced  jets  perform  well  in  the  time-averaged  sense,  although 
sometimes  at  the  expense  of  an  increase  in  unsteadiness.  Results  are  described  in  papers  [cl2  J6]. 

The  jet  in  crossflow  has  wide  application  in  industrial  and  environmental  flows.  Dilution  jets  in 
combustors  and  cooling  jets  in  turbines  are  examples  of  jets  in  crossflow  in  industrial  flows.  In  the 
atmosphere,  flows  from  smokestacks  and  chimneys  as  well  as  accidental  ground  level  releases  can  be 
modeled  as  jets  in  crossflow.  In  some  of  these  flows,  buoyancy,  heat  release  and  compressibility  may 
play  important  roles.  However,  the  large-scale  features  of  those  flows  may  have  much  in  common 
with  a  jet  in  crossflow  in  an  incompressible  flow.  This  study  focuses  on  the  study  of  mixing  in  such 
a  fundamental  setting  through  the  tool  of  direct  simulations.  By  modulating  the  jet  inflow  into  the 
domain  at  a  variety  of  frequencies,  the  effect  of  forcing  on  mixing  is  explored  through  a  variety  of 
steady  and  time-varying  measures. 

Study  of  the  jet  in  crossflow  has  grown  in  recent  years  with  new  work  by  experimental  and  com¬ 
putational  groups  around  the  world.  Detailed  flow  visualization  studies  have  been  performed  by  1 19], 
[20]  and  [32].  The  structure  of  the  near  field  of  the  jet  in  crossflow  was  studied  by  [19],  while  [20] 
focused  on  the  origin  of  vorticity  in  the  wake  of  the  jet  in  crossflow  and  found  that  separation  events 
in  the  boundary  layer  surrounding  the  jet  led  to  the  formation  of  wake  vortices.  [32]  further  studied 
the  formation  of  vorticity  in  the  wake  and  the  patterns  of  separation  around  the  jet  as  well  as  in  the 
pipe  through  which  the  jet  flows  into  the  domain. 


Numerical  Method 


The  governing  equations  used  in  these  simulations  of  the  jet  in  crossflow  are  the  incompressible  Navier- 
Stokes  equations 


dui  du{Uj 
dt  dxj 


1  d2U{  dp 

Re  dxjdxj  di{ 


P  =  0 

uXi 

An  convection-diffusion  equation  is  used  to  model  the  evolution  of  the  scalar  field: 

ds  dujS  _  1  d^s 

dt  dxj  Pr  Re  dxjdij 


(3.17) 

(3.18) 

(3.19) 
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Both  the  Reynolds  number  of  the  channel  flow  Re/,,  based  on  the  free  stream  velocity  U0 0  and  the 
channel  half-height,  and  the  Reynolds  number  Red  based  on  the  jet  diameter  d  and  the  centerline  jet 
exit  velocity  Uj  are  3000,  while  the  velocity  ratio  r  =  Uj/Uoo  is  six  in  these  simulations.  The  Prandtl 
number  is  taken  to  be  that  of  air,  i.e.  Pr  =  0.7. 

Periodic  boundary  conditions  are  imposed  in  the  streamwise  and  spanwise  directions,  and  these  di¬ 
rections  are  treated  using  the  pseudospectral  Fourier  method.  The  wall-normal  direction  is  discretized 
using  a  staggered,  second-order  finite  difference  technique  developed  by  [8]  .  In  the  wall-parallel  di¬ 
rections,  the  computation  of  the  nonlinear  terms  is  performed  on  a  grid  which  contains  3/2  as  many 
collocation  points  as  the  computational  grid  to  eliminate  the  possibility  of  aliasing  error.  The  Navier- 
Stokes  equations  are  solved  using  a  fractional  step  technique  developed  by  [5].  The  equations  are 
advanced  in  time  by  a  hybrid  implicit-explicit  scheme  where  the  wall-normal  convective  and  diffusive 
terms  are  treated  implicitly.  The  implicit  terms  are  handled  by  a  Crank-Nicolson  scheme  while  an 
explicit,  low-storage,  five-stage,  fourth-order  Runge-Kutta  scheme  [10]  is  used  for  explicit  terms  — 
the  wall-parallel  convective  and  diffusive  terms.  The  implicit  treatment  of  the  wall-normal  convection 
eases  the  restriction  on  the  time  step  that  would  be  imposed  by  the  flow  of  the  jet  across  the  grid, 
which  is  refined  in  the  wall-normal  direction  to  resolve  the  large  velocity  gradients  close  to  the  wall. 
A  linearization  of  the  wall-normal  convective  term  in  the  wall-normal  momentum  equation  which  pre¬ 
serves  the  second-order  accuracy  of  the  time  migration  [5]  is  employed  to  avoid  the  implicit  solution 
of  a  nonlinear  equation: 


1  ( 0l4+1Uo+1  \  5«2u2+1  /~\i  a  .i2\ 


(3.20) 


The  evolution  equation  for  the  scalar  field  is  advanced  in  time  using  the  same  hybrid  implicit-explicit 
scheme.  However,  the  wall-normal  direction  is  treated  with  an  upwinding  finite  volume  scheme  for  the 
scalar,  and  filtering  of  the  high  wavenumbers  in  the  wall-parallel  direction  is  also  performed  to  avoid 
instabilities  in  the  scalar  computations. 


Computational  Setting 

Although  the  use  of  a  Fourier  method  for  the  wall-parallel  directions  provides  a  high  degree  of  resolu¬ 
tion,  the  grid  is  restricted  to  be  uniform  in  these  directions  and  cannot  be  refined  in  the  neighborhood 
of  the  jet  exit.  As  a  result,  the  resolution  in  our  simulations  is  marginal  in  the  neighborhood  of  the  jet 
exit.  The  computational  grid  is  approximately  25  x  12  x  10  d  in  the  streamwise,  wall-normal  and  span- 
wise  (x,  y  and  z)  directions,  respectively.  The  grid  spacing  in  the  streamwise  and  spanwise  directions 
is  Ax  =  A z  =  0.16dj  on  the  computational  grid,  so  that  the  jet  orifice  is  represented  by  approximately 
twelve  points.  The  grid  in  the  wall-normal  direction  is  defined  using  a  hyperbolic  tangent  stretching 
function.  The  computational  grid  is  168  x  120  x  64,  with  256  x  120  x  96  collocation  points  used  in 
the  computation  of  the  nonlinear  terms  to  eliminate  aliasing  error  in  the  wall-parallel  directions. 

The  domain  of  the  simulations  is  channel  flow,  with  periodic  boundary  conditions  in  the  streamwise 
and  spanwise  directions  and  permeable,  no-slip  walls  at  the  top  and  bottom  of  the  boundary.  The 
jet  inflow  is  specified  as  part  of  the  (possibly  time-varying)  boundary  condition  on  the  wall-normal 
velocity  field.  The  boundary  condition  on  the  scalar  is  also  specified  at  the  top  and  bottom  wall, 
allowing  the  scalar  to  be  non-zero  only  at  the  jet  exit.  The  inflow  boundary  condition  upstream  of 
the  jet  is  enforced  (weakly)  through  the  use  of  a  fringe  region  which  removes  disturbances  15  to  20 
diameters  downstream  of  the  jet,  creating  a  laminar  inlet  flow  at  the  entrance  to  the  box.  The  inlet 
flow  is  laminar,  but  the  mean  velocity  profile  is  fuller  than  the  parabolic  steady  state  profile.  The 
simulations  used  a  profile  of  U(y)  =  [/«,  (1  —  y8),  which  gives  an  initial  boundary  layer  thickness  of 
(J99%  =  0.44/i  «  2d.  Thus,  the  boundary  layer  will  grow  as  the  flow  moves  through  the  domain  but 
will  not  thicken  substantially  before  reaching  the  jet  orifice. 
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The  jet  is  placed  near  the  entrance  to  the  domain  in  the  streamwise  direction  and  is  centered  in 
the  box  in  the  spanwise  direction.  The  velocity  profile  at  the  jet  exit  is  axisymmetric  and  is  specified 
as  a  function  of  the  radial  distance  from  the  centerline  of  the  jet.  A  raised  cosine  filter  is  applied  in 
Fourier  space  to  the  velocity  profile  at  the  jet  exit  to  smooth  the  profile,  preventing  the  excitation  of 
high  wavenumber  oscillations  in  the  velocity  field  by  the  jet.  A  region  of  distributed  suction  is  placed 
near  the  exit  of  the  domain  to  ensure  that  there  is  no  net  mass  flux  across  each  of  wall. 

The  jet-in-crossflow  has  been  simulated  in  a  periodic  channel-flow  geometry,  with  the  outflow 
velocity  profile  of  the  jet  specified  and  the  jet  disturbances  at  the  outflow  of  the  domain  removed  by 
a  fringe  region  to  allow  a  nearly  laminar  inflow  into  the  domain.  The  flow  in  the  pipe  leading  to  the 
jet  outlet,  which  may  separated  before  entering  the  domain  at  low  jet  vuocity  ratios  [32],  has  been 
neglected  in  these  simulations.  The  velocity  profile  where  the  pipe  exhaust',  into  the  domain  has  been 
specified  in  the  simulations  as  a  Dirichlet  boundary  condition  on  the  normal  velocity,  and  the  wall- 
parallel  velocities  on  this  plane  have  been  set  to  zero.  (The  wall-parallel  velocities  could  be  non-zero, 
for  example  in  the  case  of  a  swirling  jet-in-crossflow.)  Imposing  the  velocity  profile  does  neglect  the 
development  of  the  pipe  flow  approaching  the  jet  exhaust,  but  may  be  a  reasonable  approximation  at 
higher  jet  velocity  ratios.  [32]  noted  that  the  flow  in  the  pipe  through  which  the  jet  flows  into  the 
domain  did  not  separate  for  velocity  ratios  greater  than  six.  The  velocity  ratio  Uj/Uoo  is  six.  Although 
a  top  hat  profile  is  desirable  for  the  jet  inflow,  the  non-dissipative  Fourier  spectral  method  used  in  the 
discretization  of  the  wall-parallel  directions  requires  smoothness  in  the  velocity  boundary  conditions. 
As  a  result,  the  jet  inflow  is  specified  as  an  exponential  function  of  the  radial  distance  from  the  center 
of  the  jet  and  is  then  filtered  using  a  raised  cosine  filter  to  remove  high  wavenumber  content.  The 
functional  form  of  the  jet  inflow  before  filtering  is: 

Uj(r)  ~  exp(—  )  (3‘21) 

where  dj  is  the  specified  diameter  of  the  jet.  Note  that  the  velocity  is  not  zero  at  the  edge  of  the  jet 
(r  =  dj/ 2);  it  does,  however,  fall  to  zero  quickly  outside  the  edge  of  the  jet.  The  mass  flux  into  the 
domain  from  the  jet  is  balanced  by  weak,  distributed  suction  near  the  end  of  the  computation  domain. 
This  suction  overlaps  with  the  fringe  region  and  ensures  that  there  is  zero  net  mass  flux  across  each 
of  the  wall-parallel  planes  in  the  channel  in  our  incompressible  flow. 

Our  computational  domain  is  periodic  in  the  streamwise  and  spanwise  directions.  If  the  fluctuations 
induced  by  the  jet-in-crossflow  are  not  removed  before  the  end  of  the  domain,  they  would  substantially 
influence  the  flow  dynamics  around  the  emerging  jet.  A  fringe  region  is  added  to  the  simulation  domain 
to  remove  the  fluctuations  induced  by  the  jet-in-crossflow  and  enforce  the  desired  inflow  conditions. 
The  inflow  is  specified  as  I7oo(y)  =  (1  —  (y/h)®)  with  the  wall-normal  and  spanwise  velocities  and  the 
scalar  set  to  zero  at  the  entrance  to  the  domain.  Here,  y  €  [~h,  h]  where  h  is  the  half-height  of  the 
channel.  The  displacement  thickness  of  the  boundary  layer  (computed  over  half  of  the  channel)  is 
somewhat  smaller  than  the  diameter  of  the  jet  exit  ( 5/dj  =  0.7),  while  the  momentum  thickness  is 
about  a  third  of  the  jet  diameter  9/dj  =  0.3).  The  Reynolds  number  of  these  simulations  based  on  jet 
diameter  and  free  stream  velocity  Re<x>  =  Uoodj/v  is  500  with  a  velocity  ratio  Uj/Uoo  of  six.  One  may 
also  define  the  Reynolds  number  based  on  the  jet  exit  velocity  and  jet  diameter  Rtj  —  Ujdj/v  which 
is  3000  in  these  simulations. 

Results 

For  a  range  of  frequencies  similar  to  but  lower  than  those  for  round  jets,  forcing  has  a  substantial 
impact  on  the  structure  and  mixing  of  the  jet.  Forcing  in  the  frequency  range  (SU  =  0.1  -  0.25)  led  to 
increased  spreading  of  the  scalar  in  the  wall-normal  direction  and  greater  penetration  into  the  crossflow. 
In  these  cases,  the  mean  scalar  profiles  are  qualitatively  different  from  the  unforced  case,  appearing 
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pinched  in  the  center  with  legs  of  fluid  stretched  below.  Fluctuations  at  the  upper  edge  of  the  jet  are 
enhanced  by  the  forcing  which  excites  the  Kelvin-Helmholtz  rollup.  As  expected  (e.g.  see  [73]),  forcing 
the  jet  at  higher  frequencies  results  in  a  jet  which  resembles  closely  the  unforced  jet.  Examination  of 
the  frequency  spectra  at  different  locations  in  the  jet  shows  that  forcing  at  relatively  high  frequencies 
Std  =  0.64  persists  only  very  close  to  the  jet  exit.  Further  from  the  jet,  no  peak  in  the  spectra  is 
visible  at  the  forcing  frequency:  apparently,  the  jet  in  crossflow  is  not  receptive  to  forcing  at  higher 
frequencies.  Forcing  at  the  lower  frequencies  Std  =  0.1-0.25  results  in  increased  energy  fax  away  from 
the  jet  exit  in  this  lower  frequencies,  providing  an  excitable  band  of  frequencies  which  are  beneficial 
for  mixing.  Figure  3.16  illustrates  the  effect  of  forcing  on  a  downstream  plane  transverse  to  the  jet. 
Forcing  can  have  a  substantial  impact  on  the  structure  of  and  mixing  by  the  jet  in  crossflow.  We  have 


Time-averaged  scalar  fields:  conventional  measures  of  mixing 
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Figure  3.16: 

simulated  the  jet  in  crossflow  to  examine  the  response  of  the  jet  to  sinusoidal  forcing  at  a  variety  of 
frequencies  in  the  range  Std  =  0.1-0.64.  The  effect  of  the  forcing  on  the  jet  is  most  substantial  in  the 
range  Std  =  0.1  -  0.25,  with  increased  spreading,  penetration  and  entrainment  at  those  frequencies. 
The  appropriate  metric  for  mixing  in  the  jets  in  crossflow  is  an  open  question.  Increased  spreading 
and  entrainment  often  comes  at  the  expense  of  increased  unsteadiness  and  intermittency  in  the  jet. 
The  acceptability  of  this  unsteadiness  is  primarily  dependent  on  the  timescale  associated  with  the 
unsteadiness.  In  a  cooling  application  for  example,  unsteadiness  at  timescales  much  faster  than  the 
characteristic  heating/cooling  timescale  of  a  material  may  not  result  in  substantial  thermal  cycling. 
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However,  slower  timescales  could  lead  to  substantial  cycling  and  fatigue.  Careful  interpretation  may 
lead  to  greater  understanding  of  the  importance  of  unsteadiness  and  the  potential  for  increased  mixing 
by  forced  jets  in  crossflow. 

3.5.8  Three-dimensional  Eulerian-Lagrangian  vortex  sheet  model  for  the  actuated 
jet  in  cross  flow 

A  novel  three-dimensional  Eulerian/Lagrangian  vortex  sheet  formulation  was  proposed  to  model  the 
actuated  jet  in  cross  flow  by  Razvan  Florea.  An  earlier  version  of  a  similar  model  was  proposed  by 
Bernd  Noack  and  Andrzej  Banaszuk. 

The  jet  is  modeled  by  a  cylindrical  vortex  sheet..  Starting  with  the  incompressible  vorticity  equa¬ 
tion,  an  integral  formulation  for  the  strength  of  the  vortex  sheet  in  a  moving  coordinate  system  is 
derived.  The  velocity  field  induced  by  the  vortex  sheet  is  given  by  the  Biot-Savart  law.  The  advantage 
of  the  present  formulation  is  that  a  fixed  number  of  states  is  required  to  describe  the  evolution  of 
the  jet  in  time.  Such  a  formulation  is  very  attractive  for  further  developing  of  control  strategies  to 
enhance  molecular  mixing. 

The  motivation  for  the  present  approach  is  to  derive  a  model  that  can  both  capture  the  mixing 
characteristics  of  a  jet  in  cross  flow  and  be  further  used  in  active  control  studies  of  the  actuated  jet. 
Some  of  the  present  models  based  on  LES/DNS  are  very  sophisticated  and  require  an  extremely  large 
number  of  states  and  long  time  simulations  with  little  physical  insight.  Other  models,  physical-based 
reduced  order  models,  while  are  computational  less  expensive,  are  either  derived  in  a  Lagrangian 
frame  with  a  variable  (non-fixed)  number  of  states  [16],  [17],  or  are  too  simplified  and  less  relevant 
from  the  physical  point  of  view  [27].  For  a  review  of  the  jet  in  crossflow  research  we  refer  to  a  paper 
by  Margason  [87].  All  these  models  prove  to  be  less  appealing  for  active  control  studies.  What  we  are 
looking  for  are  reduced  order  models  that  can  capture  the  main  aspects  of  mixing  and  at  the  same 
time  be  used  in  active  control  studies.  We  want  physically  relevant  reduced  order  models  with  a  finite 
number  of  states. 

The  present  investigation  tries  to  incorporate  in  a  rigorous  manner  the  Eulerian/Lagrangian  formu¬ 
lation  into  a  novel  three-dimensional  integral  formulation  of  the  vorticity  equation.  Another  attempt 
of  the  present  formulation  is  to  maintain  the  main  characteristics  of  a  CFD  type  discretization  for  a 
relatively  smaller  number  of  states  system. 

We  are  interested  in  modeling  a  strong  jet  issuing  from  a  pipe  into  a  crossflow  using  vortex  type 
methods.  For  a  review  of  such  methods  we  refer  to  a  paper  by  Leonard  [85].  A  schematic  representation 
of  such  a  jet  is  shown  in  Fig.  3.17.  The  quantities  Uoo  and  Vjet  are  the  velocities  of  the  undisturbed 
crossflow  and  of  the  jet,  and  the  radius  of  the  jet  is  i?jet-  We  start  with  the  assumption  that  for  a 
strong  jet  in  crossflow,  the  vorticity  is  concentrated  on  the  surface  of  the  jet  [11].  Then  the  jet  can 
be  modeled  by  a  cylindrical  vortex  sheet  governed  by  the  three-dimensional  vorticity  equation  (3.11), 
the  incompressible  continuity  equation  (3.12) 

The  closure  relation  between  the  vorticity  and  velocity  is  given  by  the  vorticity-induction  equation 
(Biot-Savart  law)  (3.15).  The  velocity  v,  at  a  point  defined  by  the  position  vector  x,  is  expressed  as  the 
sum  of  an  irrotational  component,  the  incoming  crossflow  velocity  Uqo,  and  a  solenoidal  component, 
the  vorticity-induction  component  vw 

v  =  u“+v“'  v“=i/v“xF^p‘n'(xo)  (3-22) 

Here  the  quantity  xo  is  the  position  vector  of  any  point  of  the  sheet  vorticity  defining  the  jet. 

To  discretize  Eq.  (3.13)  we  consider  an  elemental  volume  V  attached  to  the  jet  moving  surface, 
and  can  move  within  the  surface.  Hence,  the  elemental  volume  V  is  not  a  material  volume,  it  has  a 
velocity  u  which  can  be  different  from  the  velocity  v  of  the  material  volume  that  coincide  with  volume 
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V  at  time  t .  The  velocity  u  is  going  to  be  the  projection  of  v  on  some  particular  local  surface  that  is 
going  to  be  defined  later.  We  integrate  Eq.  (3.13)  over  an  elemental  volume  V  and  after  applying  the 
Reynolds  transport  theorem,  we  obtain 

4  [  udV+  f  F adS  =  0  (3.23) 

at  Jv 


wdV  +  /  F 3dS  —  0 


The  quantity  S  is  the  moving  surface  describing  the  boundary  of  the  volume  V,  and  Fs  is  the  flax 
matrix  corresponding  to  the  moving  surface 

F*  =  cj(v  -  u)T  -  vuF  =  0  (3-24) 

We  assume  that,  at  any  moment  t,  the  jet  surface  can  be  parameterized  by  a  set  of  quasi-orthogonal 
system  of  coordinates  (#2,  £3),  wbere  is  along  the  jet  and  X3  is  across  the  jet  surface.  Note,  that  this 
is  a  moving  mesh  which  deforms  in  both  space  and  time  and  has  a  fixed  number  of  points.  In  addition, 
we  introduce  the  coordinate  x\  normal  to  the  jet  surface,  i.e.  form  a  three-dimensional 

system  of  coordinates,  assumed  to  be  orthogonal.  Shown  in  Fig.  3.18  are  the  mesh  parameterization 
of  the  jet  surface,  a  view  of  of  the  elemental  volume  V  =  Axi  x  ( EFHG )  =  Axi  x  Ax2  x  A23,  and 
the  local  system  of  coordinates. 

For  each  point  on  the  surface  mesh  we  denote  the  corresponding  versors  associated  with  the  local 
system  by  ri,r2  and  r3,  and  the  components  of  the  velocity  vector  by  vi,v2  and  t>3  along  these 
directions.  Note  that  the  vorticity  cj,  at  any  point  on  the  mesh,  is  tangent  to  the  surface  mesh  at  that 
point.  Hence,  u  is  defined  locally  only  by  the  two  components  u>2  and  CJ3- 

The  position  of  each  point  on  the  mesh  that  describes  the  jet  is  given  by  the  ordinary  differential 

equation 

^  =  u  =  v\T\  +  V3T3  (3.25) 

at 


Hence,  its  section  is  aloud  to  move  to  move  with  the  jet  surface  but  not  along  the  jet.  As  a  result,  the 
number  of  sections  and  states  needed  to  represent  the  jet  is  constant.  This  characteristics  defines  the 
“Eulerian”  part  of  our  Eulerian-Lagrangian  model.  In  a  three-dimensional  Lagrangian  vortex  model 
[17],  at  each  moment  in  time,  as  the  jet  penetrates  the  cross  flow,  new  solutions  and  corresponding 
states  are  introduced.  These  states  follow  the  flow  and  leave  the  domain  as  the  flow  crosses  the 
downstream  far-field  boundary.  This  transport  of  states  through  the  domain  makes  the  Lagrangian 
vortex  model  to  be  less  appealing  for  an  efficient  analysis  of  the  jet  actuation. 

Next  we  describe  shortly  the  numerical  implementation  of  our  model.  For  convenience  we  introduce 
two  quantities  72  and  73  defined  by 

J  uidV  »  ojAxi  AX2AX3  =  u>2 Aii  A1 3. A12 T2  +  u>3 A gj A13T3  (3.26) 

V  72  73 

Then  the  differential  form  of  the  Eqs.  (3.23)  and  (3.24)  in  the  local  system  of  coordinates  (xj,  12,2:3) 
is  given  by 


J^(72Ax2)  -  (Ax3)^(«273)  =  0  (3-27) 

J^(73Ax3)  +  (72Ax2)|^  +  (Ax2)^  (V273^f  "  U372)  "  (Ai3)^(u373)  =  0  (3.28) 

These  equations  are  discretized  in  space  using  centered  differences  in  the  13  direction  and  upwinded  d- 
ifferences  in  the  X2  direction.  We  integrate  numerically  Eq.  (3.15)  with  the  Gauss-Legendre  quadrature 
rule.  We  use  three-dimensional  local  cubic  interpolation  for  the  vector  x  that  describes  the  jet  surface 
and  linear  interpolation  for  the  components  72  and  73  of  u.  We  use  a  simple  corrected  Euler  (second 
order)  for  the  time  integration  Equation  (3.13)  and  the  resulting  space  discretization  of  Eqs.  (3.27) 
and  (3.28). 
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Figure  3.17:  Jet  in  cross  flow. 


We  impose  that  the  flow  around  the  wall  is  tangential  to  the  wall.  At  the  wall  surface  level,  near 
the  duct,  outside  the  jet,  the  flow  moves  around  the  jet,  tangential  to  the  jet.  We  also  impose  the 
velocity  of  the  jet  flow,  in  a  cross  section  of  the  duct,  normal  to  the  wall.  We  modeled  these  conditions 
using  rings  of  distributed  vorticity  on  the  wall,  around  the  duct,  7W3,  and  on  the  duct  wall,  bellow 
the  wall,  702  and  703.  Shown  in  Fig.  3.19  is  the  type  of  discretization  we  use  for  the  initial  jet,  duct 
and  wall.  The  above  boundary  conditions  are  solved  implicitly  at  each  moment  in  time.  Using  the 
velocity  induced  by  different  vortex  sheets  (wall,  duct  and  jet),  an  algebraic  system  of  equations  is 
formed  and  the  solution  is  found  using  an  LU  decomposition.  Note  that  the  matrix  of  this  system 
does  not  change  in  time,  the  elements  of  this  matrix  are  defined  by  the  Biot-Savart  law  and  as  a  result 
depend  only  the  geometrical  position  of  the  wall  and  duct. 

At  t  =  0,  we  impose  the  additional  condition  that  the  vorticity  distributions  below  the  wall,  around 
the  duct,  and  above  the  waU,  on  the  jet  are  identical.  This  initial  solution  is  further  marched  in  time. 

To  show  the  effectiveness  of  the  near-field  boundary  conditions,  we  show  in  Fig.  3.20  the  initial 
velocity  profiles  in  a  vertical  plane.  The  first  case  corresponds  to  a  core  size  of  0.1  of  the  local  mesh 
size,  the  second  one  corresponds  to  a  core  twice  the  local  mesh  size. 

Computing  the  correct  velocity  of  the  jet  surface  in  the  Eulerian-Lagrangian  formulation  requires 
additional  attention.  At  each  moment  in  time,  for  a  given  geometry  we  compute  the  velocity  induced 
(Biot-Savart)  by  the  vortex  sheet.  Note  that  the  surface  is  defined  by  rings  that  are  allowed  to  deform 
with  the  jet  surface  but  not  allowed  to  move  along  the  jet.  Next,  we  allow  the  vortex  surface  to  move 
in  space  and  a  new  surface  is  defined.  The  new  surface  is  interpolated  at  cross  sections  corresponding 
to  the  old  rings  and  new  rings  are  defined.  The  difference  in  space  between  the  old  and  new  rings 
determines  the  correct  velocity  of  the  jet  surface.  This  approach  requires  further  investigations  to 
understand  the  effect  of  the  additional  errors  introduced  with  the  interpolation.  These  errors  may 
become  important  when  geometrical  instabilities  form  on  the  surface  of  the  jet. 

The  finite  number  of  j  sections  which  define  the  jet,  and  the  finite  length  of  the  jet  cause  additional 
difficulties.  For  the  moment,  we  used  an  increased  core  for  the  vorticity  at  the  end  of  the  jet.  However 
the  effect  of  these  approximations  require  further  investigations. 

To  assess  the  effectiveness  of  the  near-field  boundary  conditions  we  model  the  evolution  of  a  jet  in 
cross  flow  defined  by  the  following  parameters:  Uoo  =  1.,  Vjet  =  6.,  ,  JZjet  =  1-  The  jet  is  modeled  by 
j  =  35  rings  with  k  =  19  points  per  ring. 


Figure  3.18:  Jet  in  cross  flow.  Left:  quasi-orthogonal  mesh  parameterization  of  the  jet  surface.  Right: 
elemental  volume  V  =  Aii  x  ( EFHG ),  and  local  system  of  coordinates. 


Figure  3.19:  Initial  jet,  duct  and  wall  discretization. 
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Figure  3.20:  Initial  velocity  profiles  in  a  vertical  plane  through  the  jet  center.  Left,  core  size  of  0.1  x  As. 
Right,  core  size  of  2.  x  As.  As  -  local  mesh  size  parameter. 

To  accelerate  the  convergence  of  the  jet  to  its  final  position,  the  jet  is  also  constantly  moved 
(rotated)  in  the  direction  of  the  U0 0  for  the  first  100  iterations.  We  show  in  Fig.  3.21  the  jet  evolution 
after  20,  60,  100,  200  and  1000  iterations.  In  addition,  we  also  show  in  Fig.  3.22  shape  sections  of  the 
jet  (after  1000  iterations).  We  note  that  due  to  errors  in  interpolation  (described  in  Section  C.),  the 
jet  starts  to  loose  its  symmetry. 

For  the  same  kinematic  conditions,  we  varied  initial  conditions  (length  of  jet,  imposed  motion  of 
the  jet  in  the  x  direction)  and  the  mesh  density.  We  found  that  the  most  dynamic  part  of  the  jet  is 
near  the  wall  where  the  counter  rotating  vortex  pairs  are  formed.  The  upper  part  of  the  jet  behaves  in 
a  more  static  way.  Based  on  its  length  and  inclination,  it  induces  a  component  of  the  jet  velocity  near 
the  wall  that  has  direct  effect  on  the  jet  bending  and  the  CVP  formation.  As  a  result,  the  combination 
of  the  initial  conditions  and  the  errors  in  the  interpolation  process  has  direct  effect  on  the  solution 
evolution,  preventing  the  jet  to  evolve  to  the  correct  position. 

3.F.9  Two-vortex  model  for  the  actuated  jet  in  cross  flow 

A  parabolized  two-dimensional  two-vortex  model  has  been  developed  by  Thomas  John,  Satish  Narayanan, 
Andrzej  Banaszuk,  and  Shubhro  Ghosh. 

The  dynamically  significant  flow  structures  appear  to  be  (  see  Figure  3.13)  the  shear  layer  struc¬ 
tures  generated  by  the  Kelvin-Helmholtz  (KH)  instability,  the  counter-rotating  vortex  pair  (CVP),  a 
typically  weak  horseshoe  vortex  (appearing  for  low  jet  Reynolds  numbers  near  the  jet  and  cross-flow 
bottom  juncture),  and  an  upright/wake  vortex  system  that  develops  very  close  behind  the  jet  (resem¬ 
bling  a  wake  of  vortices).  Of  these,  we  only  discuss  and  model  the  KH  and  CVP  dynamics  and  mixing 
characteristics.  In  doing  so,  we  use  inviscid  descriptions  of  the  flow  features  and  structures,  ignoring 
fine-scale  turbulence  driven  viscous  effects.  We  do  not  attempt  to  model  or  explain  the  appearance  of 
the  KH  or  CVP  structures,  but,  given  their  presence  and  from  their  observed  behavior,  we  provide  a 
modeling  framework  for  studying  the  KH-CVP  system  dynamics  to  analyze/derive  control  solutions 
for  mixing  enhancement. 

Kelso,  Lim  &  Perry  [32]  based  on  their  experimental  study  of  jets  in  cross-flow,  hypothesize  the 
formation  of  vortex  rings  (as  a  result  of  the  Kelvin-Helmholtz  instability  of  the  jet  shear  layer)  from 
the  jet  which  get  distorted  due  to  the  crossflow  such  that:  the  upstream  portion  stays  in  a  plane 
transverse  to  the  jet  and  the  downstream  portion  bends  along  the  direction  of  the  jet.  while  being 
positioned  in  a  high  speed  region  of  the  jet,  resulting  in  its  stretching  (i.e.  vorticity  intensification). 
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Figure  3.21:- Jet  evolution  after  20,  60,  100,  200  and  1000  iterations. 
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Figure  3.22:  Shape  sections  of  the  jet  at  J  =  5, 10, 15, 20, 30. 

Cortelezzi  &  Karagozian  note  similar  behavior  of  the  vortex  rings  in  their  simulations  of  a  jet  in 
cross-flow  using  3D  vortex  filament  simulations. 

Using  a  direct  numerical  simulation  (DNS)  of  jet  in  cross-flow  (P.  Blossey  &  T.Bewley,  UCSD), 
the  circulation  in  a  plane  normal  to  the  jet  trajectory  was  computed  to  obtain  the  circulation  in  the 
CVP  at  several  points  along  the  trajectory.  The  circulation  is  seen  to  increase  and  then  decrease  along 
the  trajectory  (see  Fig.  3.23). 

We  now  focus  on  the  evolution  region  from  2  to  6  jet  diameters,  where  the  KH  formation  has  been 
initiated  and  the  out-of-plane  vorticity  necessary  to  initiate  the  CVP  motion  has  already  been  set  up. 

We  simplify  and  represent  the  array  of  vortex  rings  by  “horizontal”  (transverse  to  the  jet  direction) 
rods  of  vorticity,  representing  the  upstream  KH  ring  and  two  growing  (stretching)  “legs”  to  represent 
the  jet  portion  that  has  been  bent  along  the  jet  direction(see  fig.  3.24). 

(More  features  may  be  added  if  they  are  felt  to  be  important  for  mixing  studies.)  The  transvjrse 
rod  is  moving  at  a  speed  of  U/2  (typical  convection  velocity  of  structures  in  a  shear  layer)  axJ  the 
legs  are  growing  with  respect  to  the  rod  at  a  rate  of  U/2  (since  they  are  in  a  portion  of  the  jet  which 
is  moving  at  a  higher  velocity,  say  U) 

As  mentioned  before,  the  legs  grow  away  from  the  transverse  rod  at  the  same  rate  that  the  rod 
moves  away  from  the  jet  exit.  So,  if  the  rod  if  si  away  from  the  pipe,  the  tip  of  the  leg  is  2si  from 
the  mouth  (see  fig.  3.25). 

The  long  bar  on  the  side  represents  the  initial  circulation  in  the  CVP  (To).  As  seen  from  the  figure, 
if  we  were  to  add  up  the  circulation  in  a  plane  transverse  to  the  jet,  its  value  would  increase  as  one 
moves  downstream. 

Now  we  consider  as  our  model  plane,  a  plane  moving  at  a  speed  of  U  along  the  jet  trajectory.  This 
plane  sees  an  array  of  vortex  rings  (in  reality  from  KH  structures  formed)  that  are  moving  at  a  speed 
of  U/2  into  it  (see  fig.  3.5.9). 

As  each  of  these  rods  crosses  the  plane,  it  leaves  its  leg  in  the  plane  forever  and  thus  the  circulation 
of  the  CVP  on  the  plane  is  always  increasing  ;  note  the  similar  circulation  growth  feature  from  the 
DNS  studies  summarized  earlier.  Furthermore,  note  that  the  frequency  of  KH  rods  is  twice  that  of 
the  pulsations  of  the  CVP.  Please  see  figures  3.5.9,  3.5.9  &  3.5.9. 

The  slope  of  the  increase  of  the  CVP  is  thus  dependent  on  two  quantities  -  the  frequency  of  the 
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Figure  3.29: 


KH  roll-up  and  its  strength.  Below,  we  evaluate  the  effect  of  changing  these  parameters. 

If  we  fix  the  momentum  ratio  (  jet  to  cross-stream),  the  control  parameters  are  A  which  is  the 
spacing  of  the  KH  vortices,  <f>  which  indicates  the  phase  difference  between  the  KH  crossings  and  the 
pulsation  of  the  CVP  at  a  fixed  (control)  frequency,  T0  which  denotes  the  initial  CVP  circulation 
and  <jqkh  which  is  the  initial  core  radius  for  the  KH  roll-up  (a  function  of  the  initial  boundary  layer 
thickness).  The  expressions  for  the  velocity  fields  can  be  written  as 


vx{x,y,t)  = 


SL(1_jfe)  +  Eg£: 


_.2  v 

e  cvp  ) 


(3.29) 


where 


+4?  W  (3.30) 

rl„  =  (l  -  if  +  y2  (3-31) 

&cvp  —  4*  4 vt  (3.32) 

rW-r.  +  (^+^^+W  (3-33) 

<f>  €  [0, 2ir)  (3.34) 


vy{x,y,t) 

where 


(3.35) 


r2  =  (y-/)2  +  z2  (3.36) 

Zi  =  (i  —  1)A kh  -  ^ t  (3.37) 

of  =  cr  \  Ho  +  4vt  (3.38) 

6  =  0.0025s3  -0.0453s2  +  0.1114s  +  1.4518;  s  =  Ut  (3.39) 
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Which  may  in  turn  be  represented  in  terms  of  the  control  parameters  as 

&  =  y>  t)  =  Ti(p\iP2iP3iP4ii)F(x,y,t)  (3.40) 

y  =  vy(x,y,t)  =  r){p\,p2,P3,P\,t)G{x,y,t)  +  +  Asin(0(t))  (3.41) 

where 

Pi  =  To 
P2  =  A kh 

P3  =  <t> 

P4  =  <70  KH 

F,G  &  H  are  time-dependent  non-linear  functions 
rj  contains  a  time-independent  term,  a  linear  growth  in  time  and  am  time-oscillating  term 

All  vortex  cores  were  assumed  to  be  distributed  as  Gaussian  functions,  taking  into  account  that 
their  core  is  growing  in  time  (thereby  mimicking  viscous  diffusion  and  avoiding  singularities).  A  series 
r,i '  these  KH  rods  was  first  arranged  parallel  to  each  other  in  a  plane  perpendicular  to  the  model  plane 
ind  allowed  to  cross  the  plane  with  a  velocity  of  U/2.  The  CVP  in  the  plane  is  growing  and  pulsing 
as  described  above. 

Passive  scalar  transport  is  described  by  the  advection-diffusion  equation  (3.16).  For  the  flows 
of  interest  the  advection  term  v(x,  y,t)TVc(x,y,f)  dominates  the  diffusion  term  pAc(x,y,t).  The 
advection  diffusion  equation  can  be  studied  in  an  Eulerian  frame.  Various  control  problems  (including 
optimal  ones)  will  be  studied  with  this  description,  as  various  mixing  measures  can  be  determined  by 
penalizing  nonuniformity  of  the  passive  scalar  concentration  c(x,y,t)  in  space  and  time.  In  fact,  in 
a  parabolized  model,  mixing  performance  can  be  specified  in  terms  of  of  c(x,  y,  Tfinai),  where  Tjinai 
is  the  time  when  the  transverse  plane  reaches  a  fixed  position  downstream  (representing  position  of  a 
turbine  inlet  guide  vane).  Note  that  the  result  of  the  mixing  performance  evaluation  will  depend  on 
the  initial  condition.  In  case  of  periodic  bulk  jet  forcing,  one  should  calculate  a  cumulative  measure 
of  mixing  performance  over  one  period  of  forcing  from  “instantaneous”  measures  of  performance  that 
depend  on  the  initial  conditions  corresponding  to  a  particular  control  phase. 

Alternatively,  neglecting  the  diffusion  term,  one  can  study  the  passive  scalar  mixing  using  La- 
grangian  simulation  advecting  tracer  particles  along  the  direction  of  the  flow  given  by  velocity  filed 
v(x,  y,t).  The  mixing  can  be  studied  using  methods  of  dynamical  system  theory. 

In  a  preliminary  study,  particles  were  introduced  on  the  plane  with  a  Gaussian  probability  dis¬ 
tribution  centered  at  the  origin  and  the  effect  of  KH  frequency  on  the  final  particle  distribution  was 
examined.  One  can  see  from  Figures  3.30  and  that  low  frequency  forcing  was  much  more  effective  in 
spreading  the  jet  particles  than  high  frequency.  Note  that  the  parabolized  two-vortex  model  does  not 
explain  the  creation  or  interactions  of  the  KH  or  CVP  structures,  hence,  it  does  not  satisfy  Property 
1  (physical  relevance).  However,  one  could  imagine  construction  of  a  “grey  box”  model  that  would 
use  parametrization  of  the  vorticity  field  with  temporal  coefficients  dynamics  reconstructed  from  DNS 
simulations.  The  parabolized  model  is  not  arbitrarily  refinable,  but  one  could  develop  a  refinable  3D 
version.  However,  even  in  its  present  version  the  model  is  still  very  useful,  as  it  provides  a  framework 
for  studying  effects  of  KH-CVP  structures  on  particle  mixing.  Note  that  the  model  is  suitable  for 
control  design  using  methods  of  dynamical  systems  and  control  theory.  In  particular,  one  can  study 
desirable  dynamics  of  coherent  structures  using  dynamical  systems  or  control  theory  and  limits  of 
achievable  mixing  performance.  Once  the  desired  behavior  has  been  identified,  one  can  attempt  to 
enforce  this  behavior  in  a  flow  model  (Eulerian-Lagrangian  or  DNS)  by  solving  numerically  a  tracking 
problem,  which  in  turn  can  be  approached  as  an  optimization  problem. 


(3.42) 

(3.43) 

(3.44) 

(3.45) 
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Figure  3.30:  A  snapshot  of  a  typical  velocity  field,  and  particle  trajectories  for  high  (left)  and  low 
(right)  control  frequency 


Two-Vortex  model  allows  to  study  mixing  performance 
measures  and  beneficial  flow  structures 
Results  ofLagrangian  particle  simulations  (for  zero  diffusion) 
show  that  low  frequency  is  effective  for  mixing  Rnal  distribution 


Initial  distribution  of  particles  High  frequency  jet  pulsing  Low  frequency  jet  pulsing 


2sx  2s bins,  s=l,...,7. 

Concentration  in  bin  (k,l)  at  scale  s: 


Cs  = 


nA 


u 


nAu  +  nBu 


Mixing  Variance  Coefficient  (I.  Mezic)  at  scale  s: 

MVC‘  =  2^r^ 


MVC  for  high  frequency  MVC  for  low  frequency 


Figure  3.31:  Effect  of  high  and  low  control  frequency  on  the  Mixing  Variance  Coefficient 
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3.5.10  Control  of  jet  in  cross  flow:  experimental  results 

Experiments  with  control  of  jet  in  cross  flow  were  performed  in  a  UTRC  facility  shown  in  Figure  3.32 
by  Satish  Narayanan  and  Prabir  Barooah.  This  work  was  funded  from  UTRC  internal  DCE  funds,  but 
it  was  closely  tied  to  AFOSR  effort.  Some  experiment  were  guided  by  prediction  of  AFOSR-funded 
project  described  in  this  report.  The  data  collected  from  the  experiment  will  be  used  to  validate 
models  developed  in  AFOSR  project. 


Experimental  facility  for  jet  in  cross-flow  control 


Figure  3.32:  Experimentally  measured  mean  velocity  profiles  in  jet  in  cross  flow  experiments  at  UTRC 

Objectives  of  the  experiments  with  control  of  jet  in  cross  flow  was  to  determine  physical  mechanisms 
to  be  affected  by  actuation  to  achieve  mixing  enhancement  and  demonstrate  mixing  enhancement  The 
experimental  rig  had  a  non-reacting,  single  jet  in  cross-flow.  Diagnostics  involved  hot-film  anemomtry, 
and  Mie  scattering.  High  bandwidth  actuator  was  used  to  excite  the  whole  jet  flow. 

Cold,  non- reacting  sub-scale  experiments  were  conducted  of  the  forced  and  unforced  jet  in  crossflow 
with  a  velocity  ratio  of  6  (as  in  the  DNS  simulations).  Slightly  higher  Reynolds  numbers  based  on 
channel  height  were  used  but  similar  jet  Reynolds  numbers  were  used. 

A  mechanical  actuator  (namely,  a  spinning  valve)  was  used  to  modulate  the  jet  at  several  fre¬ 
quencies  around  the  forcing  frequency  predicted  by  the  simulations  as  being  the  effective  one  (namely, 
Std  =  0.1  —  0.25).  Similar  mixing  and  entrainment/mixing  enhancements  were  noted  with  flow  visual¬ 
ization  and  quantitative  measurements  of  the  velocity  field,  confirming  and  supporting  the  effectiveness 
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of  ”  low  frequency”  actuation. 

Preliminary  results  shown  in  Figure  3.33  indicate  mixing  enhancement  at  low  forcing  frequencies, 
in  good  agreement  with  predictions  of  analysis  using  DNS  and  Two- Vortex  model  presented  in  the 
previous  sections. 


Preliminary  exploration  of  low-frequency  forcing  concept: 

entrainment  &  unsteady  mixing  enhancement 


Figure  3.33:  Experimentally  measured  mean  velocity  profiles  in  jet  in  cross  flow  experiments  at  UTRC 
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Combustion  Instability” ,  submitted  to  Automatica. 

[j5]  G.  Thdmor,  A.  Banaszuk,  “Observer-based  control  of  vortex  motion  in  a  recirculation  zone”, 
submitted  to  IEEE  Transaction  on  Controls  Systems  Technology. 

[j6]  P.  Blossey,  S.  Narayanan,  T.R.  Bewley,  “Dynamics  &  control  of  a  jet  in  cross-flow:  direct  numerical 
simulations  &  experiments”,  in  preparation. 

[j7]  P.  Mehta,  A.  Banaszuk,  “On  the  identification  of  heat  release  nonlinearity  in  the  thermoacoustic 
feedback  loop” ,  in  preparation. 


5.2  Conference  papers 

[cl]  G.S.  Copeland,  I.G.  Kevrekidis,  R.  Rico-Martinez,  Adaptive  detection  of  instabilities  and  nonlinear 
analysis  of  a  reduced-order  model  for  flutter  and  rotating  stall  in  turbomachinery,  1999  CCA,  Hawaii. 
[c2]  S.  Narayanan,  A.I.  Khibnik,  C.A.  Jacobson,  I.G.  Kevrekidis,  R.  Rico-Martinez,  and  K.  Lust.  Low 
dimensional  models  for  active  control  of  flow  separation.  1999  CCA,  Hawaii. 

[c3]  A.  Banaszuk  and  J.  Hauser,  On  control  of  planar  periodic  orbits,  1999  CDC,  December  1999, 
Phoenix. 

[C4]  R.  Smith,  A.  Banaszuk,  and  G.  Dullerud,  Model  validation  approaches  for  nonlinear  feedback 
systems  using  frequency  response  measurements,  1999  CDC,  December  1999,  Phoenix. 

[c5]  A.  Banaszuk,  Y.  Zhang,  and  C.A.  Jacobson.  Active  Control  of  Combustion  Instabilities  in  Gas 
Turbine  Engines  for  Low  Emissions.  Part  II:  Adaptive  Control  Algorithm  Development,  Demon¬ 
stration  and  Performance  Limitations,  AVT  Symposium  on  Active  Control  Technology  for  Enhanced 
Performance  Operation  Capabilities  of  Military  Aircraft,  Land  Vehicles,  and  Sea  Vehicles,  May  2000, 
Braunschweig. 
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[c6l  B.  Coller,  B.  Noack,  S.  Narayanan,  A.  Banaszuk,  and  A.  Khibnik,  “Reduced-Basis  Model  for 
Active  Separation  Control  in  a  Planar  Diffuser  Flow”,  paper  AIAA-2000-2562,  Fluids  2000,  June 
2000,  Denver. 

[c7]  A.  Banaszuk,  Y.  Zhang,  and  C.A.  Jacobson.  Adaptive  Control  of  Combustion  Instability  Using 
Extremum-Seeking,  2000  ACC,  June  2000,  Chicago. 

[c8]  I.  Mezic  and  A.  Banaszuk,  “Comparison  of  Systems  with  Complex  Behavior:  Spectral  Methods”, 
2000  CDC,  Sydney,  December  2000. 

[c9]  B.R.  Noack,  I.  Mezic,  and  A.  Banaszuk,  “Controlling  vortex  motion  and  chaotic  advection”,  2000 
CDC,  Sydney,  December  2000. 

[clO]  K.B.  Ariyur,  A.  Banaszuk,  M.  Krstic,  “Identification  of  averaged  dynamics  of  a  controlled  com¬ 
bustion  instability”,  2000  CDC,  Sydney,  December  2000. 

[di]  G.  Tadmor,  A.  Banaszuk,  “Observer-based  control  of  vortex  motion  in  a  recirculation  zone”, 
submitted  to  2001  NOLCOS,  St.  Petersburg,  July  2001. 

[cl2]  P.  Blossey,  S.  Narayanan,  T.R.  Bewley,  “Dynamics  &  control  of  a  jet  in  cross-flow:  direct 
numerical  simulations  &  experiments”,  submitted  to  2001  (IUTAM)  Symposium  on  Turbulent  Mixing 
and  Combustion,  June  3-6,  2001,  Kingston. 

[C13]  p.  Mehta,  A.  Banaszuk,  “On  the  identification  of  heat  release  nonlinearity  in  the  thermoacoustic 
feedback  loop”,  submitted  to  2001  CDC,  Orlando  December  2001. 

5.3  Invited  Sessions 

The  following  invited  sessions  were  organized  with  AFOSR  support  and  contain  AFOSR-funded  pa- 
pers: 

Control  applications  in  flows  and  turbomachines  (organized  by  Scott  Copeland  and  Satish  Narayanan). 
CCA  99,  Hawaii. 

Theory  and  Applications  of  Extremum-Seeking  Control  (organized  by  Andrzej  Banaszuk  and  Mario 
Rotea).  ACC  2000,  Chicago. 

Control  of  Mixing  in  Shear  Flows  (organized  by  Andrzej  Banaszuk  and  Luca  Cortelezzi) .  CDC  2000, 
Sydney. 

Nonlinear  Model  Validation  (organized  by  Roy  Smith  and  Andrzej  Banaszuk).  CDC  2000,  Sydney. 


5.4  Talks  at  meetings 

The  talks  at  meetings  that  do  not  have  proceedings  include: 

Talk  at  Flow  Control  Worskshop,  Boston  University,  October  1999 

Andrzej  Banaszuk,  Bernd  R.  Noack,  Igor  Mezic 
A  hierarchy  of  vortex-based  models  for  coherent-structure  control 

The  shear  flow  phenomena  strongly  affect  operations  of  industrial  products  like  jet  engines,  helicopter- 
s,  and  HVAC  devices.  For  instance,  in  the  case  of  airfoils  at  high  angle  of  attack  the  flow  separates 
close  to  the  leading  edge  causing  decrease  in  lift  and  increase  in  drag.  Full  attachment  of  the  flow  is 
desirable,  but  not  yet  technologically  feasible  goal  in  most  industrial  cases,  because  of  limited  actuator 
authority.  However,  experiments  show  that  one  can  significantly  increase  the  lift  and  reduce  drag  by 
exciting  beneficial  coherent  structures  in  the  separated  shear  layer  using  relatively  small  actuation 
effort.  Abundance  of  shear  flow  phenomena  in  industrial  products,  high  potential  benefits  from  shear 
flow  control,  technical  feasibility  of  shear  flow  control  demonstrated  in  experiments  indicate  need  to  fo¬ 
cus  research  in  the  area  of  shear  flow  control.  In  our  talk  we  present  a  general  framework  for  modeling, 
analysis,  and  control  of  shear  flows  using  tools  from  dynamical  systems  and  control  theory.  We  show  an 
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example  of  application  of  the  framework  to  control  of  mixing  in  combustors.  We  also  identify  areas  of 
basic  research  in  control  of  shear  flows  that  can  support  technology  transition  in  reasonably  short  time. 

Talks  at  APS  (Americam  Physical  Society)  meeting,  November  2000,  Washington  DC 

Peter  Blossey,  Thomas  Bewley  (University  of  California,  San  Diego,  CA),  Satish  Narayanan  (United 

Technologies  Research  Center,  East  Hartford,  CT) 

Direct  numerical  simulations  of  a  jet  in  crossflow 

The  mixing  characteristics  of  a  jet  in  cross-flow  can  be  substantially  enhanced  using  unsteady  mod¬ 
ulation  of  the  jet  flow,  as  noted  in  prior  experimental  investigations.  We  explore  the  dynamics  of  a 
jet  in  crossflow  with  direct  numerical  simulations  to  assess  the  impact  of  unsteady  forcing  on  its  spa- 
tiotemporal  mixing  properties.  The  Reynolds  number  (Re  =  Ujd/u)  for  the  jet  is  3000  with  a  velocity 
ratio  ( Uj/Uoo )  of  6.  The  simulations  are  performed  in  a  periodic  channel  flow  geometry  and  include  a 
fringe  region  to  remove  disturbances  from  the  flow  at  the  domain  entrance.  In  addition  to  exploring 
the  vorticity  dynamics  in  the  flow  field,  we  study  the  statistical  properties  of  this  highly  3D,  unsteady 
flow,  particularly  of  its  cross-stream  mixing.  A  conserved  scalar  has  been  added  to  the  simulations  to 
enable  quantification  of  the  jet  mixing,  and  a  variety  of  mixing  metrics  are  assessed  to  evaluate  the 
effect  of  unsteady  jet  forcing.  Conventional,  time-averaged  measures  of  mixing  can  be  misleading  in 
this  flow.  Higher  order  statistics  of  mixing,  however,  are  shown  to  be  insightful  in  discriminating  the 
extent  of  mixing  among  various  unsteady  flow  patterns. 

Thomas  John  (Cornell  University,  Ithaca,  NY),  Satish  Narayanan,  Andrzej  Banaszuk,  Shubhro  Ghosh 
(United  Technologies  Research  Center,  East  Hartford,  CT) 

Reduced-order  model  for  mixing  in  a  jet  in  cross-flow 

We  present  a  reduced-order  model  describing  the  unsteady  mixing  associated  with  the  three-dimensional 
coherent  structure  dynamics  of  a  jet  in  cross-flow.  The  dynamically  significant  flow  features  are  the 
shear  layer  structures  emanating  in  the  jet  near  field,  and  a  pair  of  counter-rotating  vortices  emerging 
in  the  intermediate  and  far  fields  of  the  jet  (where  the  jet  begins  to  bend  and  align  with  the  oncoming 
cross  flow).  We  postulate  a  ’’parabolized”  vortex  model  consistent  with  the  observed  dominant  vortex 
dynamics.  The  velocity  and  vorticity  fields  obtained  from  the  model  compare  well  qualitatively  with 
prior  experimental  and  direct  numerical  simulation  results.  Lagrangian  particle  simulations  are  used 
to  quantify  the  extent  of  mixing,  relative  to  a  mixing  dynamical  system,  at  ieyeral  spatial  scales.  The 
model  appears  to  be  suitable  for  use  in  the  development  of  novel  control  la^j  tor  mixing  enhancement. 

Andrzej  Banaszuk  (United  Technologies  Research  Center),  Igor  Mezic  (Division  of  Engineering  and 
Applied  Sciences,  Harvard  University) 

Control  of  point  vortex  systems 

We  present  results  on  controllability  and  observability  of  point  vortex  motion  actuated  by  external 
potential  fields.  Sufficient  conditions  for  transformation  of  a  point  vortex  system  into  the  so-called  flat 
coordinates  which  allow  for  simple  description  of  controllability  and  observability  of  the  system  will  be 
presented.  Control  of  statistical,  large-scale  properties  of  point  vortex  systems  is  also  considered.  We 
present  a  number  of  examples  of  specific  vortex  configurations  such  as  counterrotating  vortex  pairs. 
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